From d3e8e7669eb0a4162484914497e10dbf86fda0eb Mon Sep 17 00:00:00 2001 From: Felix Moncada Date: Fri, 9 Feb 2024 18:16:00 +0100 Subject: [PATCH] Changes to perform SCF core ionized states with the goal of doing XES: Ionize selected orbitals with fractional occupation zero - to perform Corrected "ionizeSpecies" spelling in input Separated single SCF iterate into smaller subroutines Updated exchange orbitals convergence method - required to ensure ionization of the desied state Added core ionized H2O molecule example On going: NOCI Franck Condon factors and transition dipole integrals Matrix power and other functions computed from SVD Separated NOCI generate densities into smaller subroutines Small changes to print/write/read different molecular systems Integer Vector read/write to file Other: One body energy components for CI density Fixed output printing densities below 1E-100 --- bin/lowdin | 6 +- examples/FCI_Charry2018/H-.FCI-DZ.lowdin | 56 + examples/FCI_Charry2018/PsH.FCI-DZ.lowdin | 58 + examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin | 64 + .../GRIBAKIN-POTENTIAL/generatePotential.py | 8 +- src/CI/ConfigurationInteraction.f90 | 80 +- src/CalcProp/CalculateProperties.f90 | 25 +- src/DFT/DFT.f90 | 4 +- src/DFT/DensityFunctionalTheory.f90 | 1 + src/NOCI/NOCI.f90 | 12 +- src/NOCI/NonOrthogonalCI.f90 | 1414 ++++++++++++----- src/PT/PropagatorTheory.f90 | 156 +- src/core/CONTROL.f90 | 77 +- src/core/Matrix.f90 | 96 +- src/core/MolecularSystem.f90 | 444 +++--- src/core/Vector.f90 | 426 ++++- .../IntegralTransformation.f90 | 14 +- .../TransformIntegralsC.f90 | 68 +- .../TransformIntegralsE.f90 | 74 +- src/ints/Ints.f90 | 3 - src/output/InputOutput.f90 | 8 +- src/output/OutputBuilder.f90 | 101 +- src/scf/DensityMatrixSCFGuess.f90 | 24 +- src/scf/MultiSCF.f90 | 81 +- src/scf/SingleSCF.f90 | 238 ++- src/scf/WaveFunction.f90 | 68 +- test/H2O.APMO.UP2.IS.IMO1.IT-E.lowdin | 2 +- test/H2O.APMO.UP2.IS.IMO1.lowdin | 2 +- test/H2O.APMO.UP2.IS.IMO2.IT-E.lowdin | 2 +- test/H2O.APMO.UP2.IS.IMO2.lowdin | 2 +- test/H2O.APMO.UP2.IS.lowdin | 2 +- test/H2O.APMO.UP2.IS2.IT-E.lowdin | 2 +- test/H2O.APMO.UP2.IS2.lowdin | 2 +- test/H2O.BOA.UP2.IT-E.lowdin | 2 +- test/H2O.BOA.UP2.lowdin | 2 +- test/H2O.coreHole.lowdin | 29 + test/H2O.coreHole.py | 46 + test/HCN.APMO.RKS.parallel.lowdin | 35 + test/HCOOPs.HF.DensCube.lowdin | 3 +- test/HCOOPs.HF.DensCube.py | 6 +- test/LiF2-.EDA.lowdin | 21 + test/LiF2-.MP2.lowdin | 21 + test/LiF2-.addParticles.lowdin | 20 + test/LiF2-e+.EDA.lowdin | 21 + test/LiF2-e+.MP2.lowdin | 21 + test/LiF2-e+.ghost.lowdin | 21 + test/LiH.APMO.TOP2.IT-E.lowdin | 2 +- test/LiH.APMO.TOP2.lowdin | 2 +- test/LiH.TOP2.IT-E.lowdin | 2 +- test/LiH.TOP2.lowdin | 2 +- test/Mg.e+.gribakin.lowdin | 26 + test/Mg2.e+.gribakin.lowdin | 29 + test/PsF-RENP3.lowdin | 2 +- test/PsH.FCI.test.lowdin | 33 + 54 files changed, 2915 insertions(+), 1051 deletions(-) create mode 100644 examples/FCI_Charry2018/H-.FCI-DZ.lowdin create mode 100644 examples/FCI_Charry2018/PsH.FCI-DZ.lowdin create mode 100644 examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin create mode 100644 test/H2O.coreHole.lowdin create mode 100644 test/H2O.coreHole.py create mode 100644 test/HCN.APMO.RKS.parallel.lowdin create mode 100644 test/LiF2-.EDA.lowdin create mode 100644 test/LiF2-.MP2.lowdin create mode 100644 test/LiF2-.addParticles.lowdin create mode 100644 test/LiF2-e+.EDA.lowdin create mode 100644 test/LiF2-e+.MP2.lowdin create mode 100644 test/LiF2-e+.ghost.lowdin create mode 100644 test/Mg.e+.gribakin.lowdin create mode 100644 test/Mg2.e+.gribakin.lowdin create mode 100644 test/PsH.FCI.test.lowdin diff --git a/bin/lowdin b/bin/lowdin index 04be83d0..bfd955d0 100755 --- a/bin/lowdin +++ b/bin/lowdin @@ -285,7 +285,7 @@ if [ $extFile="lowdin" ]; then '{ printf("\n&Output\n") for(i=1;i<=NF;i++){ - if($i=="specie"){ + if($i=="species"){ printf("\tOutput_species = %s\n",$toupper((i+2)) ) } else if($i=="orbital"){ @@ -361,7 +361,8 @@ if [ $extFile="lowdin" ]; then cp $nameFile*.dens $LOWDIN_SCRATCH/$nameFile &> /dev/null cp $nameFile*.sup $LOWDIN_SCRATCH/$nameFile &> /dev/null cp $nameFile*.der $LOWDIN_SCRATCH/$nameFile &> /dev/null - cp $nameFile*.NOCI.coords $LOWDIN_SCRATCH/$nameFile &> /dev/null + cp $nameFile*.NOCI* $LOWDIN_SCRATCH/$nameFile &> /dev/null + cp $nameFile*.refNOCI* $LOWDIN_SCRATCH/$nameFile &> /dev/null mv $nameFile*.dens $LOWDIN_SCRATCH/$nameFile &> /dev/null mv $nameFile*.coup $LOWDIN_SCRATCH/$nameFile &> /dev/null mv $nameFile*.over $LOWDIN_SCRATCH/$nameFile &> /dev/null @@ -444,6 +445,7 @@ if [ $extFile="lowdin" ]; then mv $LOWDIN_SCRATCH/$nameFile/*.plainvec $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 mv $LOWDIN_SCRATCH/$nameFile/*.cub $currentPath &> 2 mv $LOWDIN_SCRATCH/$nameFile/*.dens $currentPath &> 2 mv $LOWDIN_SCRATCH/$nameFile/*.eps $currentPath &> 2 diff --git a/examples/FCI_Charry2018/H-.FCI-DZ.lowdin b/examples/FCI_Charry2018/H-.FCI-DZ.lowdin new file mode 100644 index 00000000..9659ad2c --- /dev/null +++ b/examples/FCI_Charry2018/H-.FCI-DZ.lowdin @@ -0,0 +1,56 @@ +!The goal of this calculation is to compute the binding energy of a positron bound complex +!Reported: +!E(H-): -0.524029 + +SYSTEM_DESCRIPTION='H- from Charry 2018 (10.1002/anie.201800914)' + +GEOMETRY + e-(H) AUG-CC-PVDZ 0.00 0.00 0.00 addParticles=1 + H dirac 0.00 0.00 0.00 +END GEOMETRY + +!method to solve the SCF - CI only works for unrestricted reference +!CI level strings chooses the desired excitations to be included. FCI is all possible excitations + +TASKS + method = "UHF" + configurationInteractionLevel ="FCI" + !configurationInteractionLevel ="CIS","CISD","CISDT","CISDTQ" +END TASKS + +!Compute only the "numberOfCIstates" states. Here we only need the ground state +!Compute the density matrix for "CIstatesToPrint" states, for density outputs +!Generate the natural orbitals, for visualization in molden files +!The Davidson diagonalization implemented in JADAMILU is the recomended method. +!For small systems, full matrix diagonalization with DSYEVX is possible +!CI EigenVectors with coefficient higher than "CIPrintThreshold" are printed +!Printing format "OCCUPIED" shows the coefficients, "ORBITALS" shows the strings, "NONE" skips printing +!Strict SCF convergence improves the quality of the CI results (not required for the FCI) + +CONTROL + numberOfCIstates=1 + CIStatesToPrint=1 + CINaturalOrbitals=T + CIdiagonalizationMethod = "JADAMILU" + !CIdiagonalizationMethod = "DSYEVX" + !CIdiagonalizationMethod = "ARPACK" + CIPrintEigenVectorsFormat = "OCCUPIED" !"NONE","ORBITALS" + CIPrintThreshold = 5e-2 + !totalEnergyTolerance=1E-12 +END CONTROL + +!INPUT_CI block help us define the frozen core and active virtuals orbitals. Here we are not restricting the excitation space +INPUT_CI + species="E-ALPHA" core=0 active=0 + species="E-BETA" core=0 active=0 +END INPUT_CI + +!With CI, moldenFiles, 1D and 2D density slices and density cubes are good ways to visualize the density results +OUTPUTS + moldenFile state=1 + densityPlot dimensions=2 point1=0.0 0.0 -6.0 point2=0.0 0.0 6.0 state=1 +END OUTPUTS + + + + diff --git a/examples/FCI_Charry2018/PsH.FCI-DZ.lowdin b/examples/FCI_Charry2018/PsH.FCI-DZ.lowdin new file mode 100644 index 00000000..7bff2499 --- /dev/null +++ b/examples/FCI_Charry2018/PsH.FCI-DZ.lowdin @@ -0,0 +1,58 @@ +!The goal of this calculation is to compute the binding energy of a positron bound complex +!Reported: +!E(PsH): -0.734559 + +SYSTEM_DESCRIPTION='PsH from Charry 2018 (10.1002/anie.201800914)' + +GEOMETRY + e-(H) AUG-CC-PVDZ 0.00 0.00 0.00 addParticles=1 + e+ E+-H-AUG-CC-PVDZ 0.00 0.00 0.00 + H dirac 0.00 0.00 0.00 +END GEOMETRY + +!method to solve the SCF - CI only works for unrestricted reference +!CI level strings chooses the desired excitations to be included. FCI is all possible excitations + +TASKS + method = "UHF" + configurationInteractionLevel ="FCI" + !configurationInteractionLevel ="CIS","CISD","CISDT","CISDTQ" +END TASKS + +!Compute only the "numberOfCIstates" states. Here we only need the ground state +!Compute the density matrix for "CIstatesToPrint" states, for density outputs +!Generate the natural orbitals, for visualization in molden files +!The Davidson diagonalization implemented in JADAMILU is the recomended method. +!For small systems, full matrix diagonalization with DSYEVX is possible +!CI EigenVectors with coefficient higher than "CIPrintThreshold" are printed +!Printing format "OCCUPIED" shows the coefficients, "ORBITALS" shows the strings, "NONE" skips printing +!Strict SCF convergence improves the quality of the CI results (not required for the FCI) + +CONTROL + numberOfCIstates=1 + CIStatesToPrint=1 + CINaturalOrbitals=T + CIdiagonalizationMethod = "JADAMILU" + !CIdiagonalizationMethod = "DSYEVX" + !CIdiagonalizationMethod = "ARPACK" + CIPrintEigenVectorsFormat = "OCCUPIED" !"NONE","ORBITALS" + CIPrintThreshold = 5e-2 + !totalEnergyTolerance=1E-12 +END CONTROL + +!INPUT_CI block help us define the frozen core and active virtuals orbitals. Here we are not restricting the excitation space +INPUT_CI + species="E-ALPHA" core=0 active=0 + species="E-BETA" core=0 active=0 + species="POSITRON" core=0 active=0 +END INPUT_CI + +!With CI, moldenFiles, 1D and 2D density slices and density cubes are good ways to visualize the density results +OUTPUTS + moldenFile state=1 + densityPlot dimensions=2 point1=0.0 0.0 -6.0 point2=0.0 0.0 6.0 state=1 +END OUTPUTS + + + + diff --git a/examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin b/examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin new file mode 100644 index 00000000..0e02adc9 --- /dev/null +++ b/examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin @@ -0,0 +1,64 @@ +!The goal of this calculation is to compute the binding energy of a positron bound complex +!Reported: +!E(e+H2^2-): -1.279680 a.u. + +SYSTEM_DESCRIPTION='PsH from Charry 2018 (10.1002/anie.201800914)' + +!add two electrons (one for each hydrogen anion) +!remove one positron +GEOMETRY + e-(H) AUG-CC-PVDZ 0.00 0.00 -1.6 addParticles=1 + e-(H) AUG-CC-PVDZ 0.00 0.00 1.6 addParticles=1 + e+ E+-H-AUG-CC-PVDZ 0.00 0.00 -1.6 + e+ E+-H-AUG-CC-PVDZ 0.00 0.00 1.6 addParticles=-1 + H dirac 0.00 0.00 -1.6 + H dirac 0.00 0.00 1.6 +END GEOMETRY + +!method to solve the SCF - CI only works for unrestricted reference +!CI level strings chooses the desired excitations to be included. FCI is all possible excitations + +TASKS + method = "UHF" + configurationInteractionLevel ="FCI" + !configurationInteractionLevel ="CIS","CISD","CISDT","CISDTQ" +END TASKS + +!Compute only the "numberOfCIstates" states. Here we only need the ground state +!Compute the density matrix for "CIstatesToPrint" states, for density outputs +!Generate the natural orbitals, for visualization in molden files +!The Davidson diagonalization implemented in JADAMILU is the recomended method. +!For small systems, full matrix diagonalization with DSYEVX is possible +!CI EigenVectors with coefficient higher than "CIPrintThreshold" are printed +!Printing format "OCCUPIED" shows the coefficients, "ORBITALS" shows the strings, "NONE" skips printing +!Strict SCF convergence improves the quality of the CI results (not required for the FCI) + +CONTROL + numberOfCIstates=1 + CIStatesToPrint=1 + CINaturalOrbitals=T + CIdiagonalizationMethod = "JADAMILU" + !CIdiagonalizationMethod = "DSYEVX" + CIPrintEigenVectorsFormat = "OCCUPIED" + !CIPrintEigenVectorsFormat = "NONE","ORBITALS" + CIPrintThreshold = 5e-2 + !totalEnergyTolerance=1E-12 +END CONTROL + +!INPUT_CI block help us define the frozen core and active virtuals orbitals. Here we are not restricting the excitation space +INPUT_CI + species="E-ALPHA" core=0 active=0 + species="E-BETA" core=0 active=0 + species="POSITRON" core=0 active=0 +END INPUT_CI + +!With CI, moldenFiles, 1D and 2D density slices and density cubes are good ways to visualize the density results +OUTPUTS + moldenFile state=1 + densityPlot dimensions=2 point1=0.0 0.0 -6.0 point2=0.0 0.0 6.0 state=1 + densityPlot dimensions=3 point1=0.0 -6.0 -6.0 point2=0.0 -6.0 6.0 point3=0.0 6.0 -6.0 state=1 +END OUTPUTS + + + + diff --git a/lib/potentials/GRIBAKIN-POTENTIAL/generatePotential.py b/lib/potentials/GRIBAKIN-POTENTIAL/generatePotential.py index e80733e6..ba4977c1 100644 --- a/lib/potentials/GRIBAKIN-POTENTIAL/generatePotential.py +++ b/lib/potentials/GRIBAKIN-POTENTIAL/generatePotential.py @@ -25,7 +25,7 @@ polarizabilities = [(angstromToBohr**3)*x for x in polarizabilities] # a.u. maxR = int(maxR/step) -print maxR +print(maxR) angularMoment = list() for i in range(0,len(exponents)) : @@ -45,14 +45,14 @@ ii = 0 for atom in range(0,len(origin)): - print origin[atom] + print(origin[atom]) realPotentialFileName = potentialName + "."+ str(origin[atom][2]) +".data" alpha = polarizabilities[atom] r = origin[atom][2] realPotentialFile = open (realPotentialFileName, "w") - for i in xrange(1,maxR+1,1): + for i in range(1,maxR+1,1): i = i*step realPotentialFile.write(str(i) + " "+ str( -1.0*alpha/(2.0*(i-r)**4) * (1 - math.exp(-((i-r)**6)/(cutoff**6))) ) + "\n" ) realPotentialFile.close() @@ -108,7 +108,7 @@ gnuplotOutput.close() if len(coeff) == 0: - print "Fitting error, results not found" + print("Fitting error, results not found") continue coefficients = [float(k) for k in coeff] diff --git a/src/CI/ConfigurationInteraction.f90 b/src/CI/ConfigurationInteraction.f90 index b37d1d97..c1cf1304 100644 --- a/src/CI/ConfigurationInteraction.f90 +++ b/src/CI/ConfigurationInteraction.f90 @@ -1142,6 +1142,7 @@ subroutine ConfigurationInteraction_densityMatrices() character(50) :: file, wfnfile, speciesName, auxstring character(50) :: arguments(2) type(matrix), allocatable :: coefficients(:), atomicDensityMatrix(:,:), ciDensityMatrix(:,:), auxDensMatrix(:,:) + type(matrix), allocatable :: kineticMatrix(:), attractionMatrix(:), externalPotMatrix(:) integer numberOfSpecies type(matrix) :: auxdensityEigenVectors @@ -1213,8 +1214,13 @@ subroutine ConfigurationInteraction_densityMatrices() write (*,*) "==============================" write (*,*) "" - allocate( coefficients(numberOfSpecies), atomicDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), & - ciDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), auxDensMatrix(numberOfSpecies,ConfigurationInteraction_instance%nproc) ) + allocate( coefficients(numberOfSpecies), & + kineticMatrix(numberOfSpecies), & + attractionMatrix(numberOfSpecies), & + externalPotMatrix(numberOfSpecies), & + atomicDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), & + ciDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), & + auxDensMatrix(numberOfSpecies,ConfigurationInteraction_instance%nproc) ) wfnFile = "lowdin.wfn" wfnUnit = 20 @@ -1228,33 +1234,45 @@ subroutine ConfigurationInteraction_densityMatrices() ! numberOfOrbitals = ConfigurationInteraction_instance%numberOfOrbitals%values(species) numberOfOccupiedOrbitals = ConfigurationInteraction_instance%numberOfOccupiedOrbitals%values(species) - arguments(2) = speciesName - arguments(1) = "COEFFICIENTS" - ! print *, "trolo", numberOfOrbitals, numberOfContractions, numberOfOccupiedOrbitals + arguments(2) = speciesName + ! print *, "trolo", numberOfOrbitals, numberOfContractions, numberOfOccupiedOrbitals - coefficients(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & - columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2)) + arguments(1) = "COEFFICIENTS" + coefficients(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & + columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2)) - ! print *, "trololo" + arguments(1) = "KINETIC" + kineticMatrix(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & + columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2)) + + arguments(1) = "ATTRACTION" + attractionMatrix(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & + columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2)) + + arguments(1) = "EXTERNAL_POTENTIAL" + if( CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & + externalPotMatrix(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & + columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2)) + ! print *, "trololo" - do state=1, CONTROL_instance%CI_STATES_TO_PRINT + do state=1, CONTROL_instance%CI_STATES_TO_PRINT - call Matrix_constructor ( ciDensityMatrix(species,state) , & - int(numberOfContractions,8), & - int(numberOfContractions,8), 0.0_8 ) + call Matrix_constructor ( ciDensityMatrix(species,state) , & + int(numberOfContractions,8), & + int(numberOfContractions,8), 0.0_8 ) - do k=1, numberOfOccupiedOrbitals - ciDensityMatrix(species,state)%values( k, k)=1.0_8 - end do + do k=1, numberOfOccupiedOrbitals + ciDensityMatrix(species,state)%values( k, k)=1.0_8 + end do - end do + end do - do n=1, ConfigurationInteraction_instance%nproc + do n=1, ConfigurationInteraction_instance%nproc - call Matrix_constructor ( auxDensMatrix(species,n) , & - int(numberOfContractions,8), & - int(numberOfContractions,8), 0.0_8 ) - end do + call Matrix_constructor ( auxDensMatrix(species,n) , & + int(numberOfContractions,8), & + int(numberOfContractions,8), 0.0_8 ) + end do end do close(wfnUnit) @@ -1566,7 +1584,25 @@ subroutine ConfigurationInteraction_densityMatrices() end do end do - + write(*,*) "" + write(*,*) "===============================" + write(*,*) " ONE BODY ENERGY CONTRIBUTIONS:" + write(*,*) "" + do state=1, CONTROL_instance%CI_STATES_TO_PRINT + write(*,*) " STATE: ", state + do species=1, molecularSystem_instance%numberOfQuantumSpecies + write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(species)%name ) // & + " Kinetic energy = ", sum(transpose(atomicDensityMatrix(species,state)%values)*kineticMatrix(species)%values) + write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(species)%name ) // & + "/Fixed interact. energy = ", sum(transpose(atomicDensityMatrix(species,state)%values)*attractionMatrix(species)%values) + if( CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & + write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(species)%name) // & + " Ext Pot energy = ", sum(transpose(atomicDensityMatrix(species,state)%values)*externalPotMatrix(species)%values) + print *, "" + end do + print *, "" + end do + !! Natural orbitals if (CONTROL_instance%CI_NATURAL_ORBITALS) then diff --git a/src/CalcProp/CalculateProperties.f90 b/src/CalcProp/CalculateProperties.f90 index cd794015..2994f1ba 100755 --- a/src/CalcProp/CalculateProperties.f90 +++ b/src/CalcProp/CalculateProperties.f90 @@ -293,8 +293,6 @@ subroutine CalculateProperties_showPopulationAnalyses(this) do speciesID = 1, MolecularSystem_getNumberOfQuantumSpecies() do type= 1, size(analysis) - - if(analysis(type) .eq. "LOWDIN" .and. CONTROL_instance%NONORTHOGONAL_CONFIGURATION_INTERACTION) cycle speciesName = trim(MolecularSystem_getNameOfSpecie( speciesID )) @@ -431,14 +429,14 @@ function CalculateProperties_getPopulation( this, typeOfPopulation, speciesID, t auxMatrix%values = matmul(this%densityMatrix(speciesID)%values, this%overlapMatrix(speciesID)%values ) - auxMatrix = Matrix_pow( this%overlapMatrix(speciesID), 0.5_8 ) + auxMatrix = Matrix_pow( this%overlapMatrix(speciesID), 0.5_8, method="SVD" ) auxMatrixB = auxMatrix auxMatrix%values = matmul( matmul( auxMatrixB%values , this%densityMatrix(speciesID)%values), auxMatrixB%values ) if(trim(speciesName) .eq. "E-ALPHA") then otherAuxMatrix%values = matmul(this%densityMatrix(otherSpeciesID)%values, this%overlapMatrix(otherSpeciesID)%values ) - otherAuxMatrix = Matrix_pow( this%overlapMatrix(otherSpeciesID), 0.5_8 ) + otherAuxMatrix = Matrix_pow( this%overlapMatrix(otherSpeciesID), 0.5_8, method="SVD" ) otherAuxMatrixB = otherAuxMatrix otherAuxMatrix%values = matmul( matmul( otherAuxMatrixB%values , this%densityMatrix(otherSpeciesID)%values), otherAuxMatrixB%values ) end if @@ -532,12 +530,31 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this) print *,"" print *," ELECTROSTATIC MOMENTS:" print *,"======================" + print *,"" + print *,"" + print *,"DIPOLE: (A.U.)" + print *,"------" + print *,"" + write (6,"(T19,4A13)") "","", ""," |D|" + + do i=1, numberOfSpecies + dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecie(this, i) + totalDipole(:)=totalDipole(:)+dipole(i,:) + write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecie( i )), dipole(i,:) + end do + dipole(numberOfSpecies+1,:)=CalculateProperties_getDipoleOfPuntualCharges() + totalDipole(:)=totalDipole(:)+dipole(numberOfSpecies+1,:) + write (6,"(T5,A15,3F13.8)") "Point charges: ", dipole(numberOfSpecies+1,:) + write (6,"(T22,A28)") "___________________________________" + write (6,"(T5,A15,3F13.8, F13.8)") "Total Dipole:", totalDipole(:), sqrt(sum(totalDipole(:)**2.0 ) ) + print *,"" print *,"DIPOLE: (DEBYE)" print *,"------" print *,"" write (6,"(T19,4A13)") "","", ""," |D|" + totalDipole=0.0_8 do i=1, numberOfSpecies dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecie(this, i)*2.54174619 totalDipole(:)=totalDipole(:)+dipole(i,:) diff --git a/src/DFT/DFT.f90 b/src/DFT/DFT.f90 index 757f46a4..218a14bc 100644 --- a/src/DFT/DFT.f90 +++ b/src/DFT/DFT.f90 @@ -68,10 +68,10 @@ program DFT select case ( job ) case ("BUILD_SCF_GRID") call DensityFunctionalTheory_buildSCFGrid() - return + STOP case ("BUILD_FINAL_GRID" ) call DensityFunctionalTheory_buildFinalGrid() - return + STOP end select !!!Computing energy and potential jobs diff --git a/src/DFT/DensityFunctionalTheory.f90 b/src/DFT/DensityFunctionalTheory.f90 index 251eea8b..97f0197b 100644 --- a/src/DFT/DensityFunctionalTheory.f90 +++ b/src/DFT/DensityFunctionalTheory.f90 @@ -88,6 +88,7 @@ subroutine DensityFunctionalTheory_buildFinalGrid() print *, "" end if call GridManager_buildGrids( "FINAL" ) + if(CONTROL_instance%GRID_STORAGE .ne. "DISK") call Functional_createFunctionals( ) if (CONTROL_instance%GRID_STORAGE .eq. "DISK") then call GridManager_writeGrids( "FINAL" ) call GridManager_atomicOrbitals( "WRITE","FINAL" ) diff --git a/src/NOCI/NOCI.f90 b/src/NOCI/NOCI.f90 index 9fa8ba3a..531caa8b 100644 --- a/src/NOCI/NOCI.f90 +++ b/src/NOCI/NOCI.f90 @@ -74,7 +74,11 @@ program NOCI if(.not. CONTROL_instance%ONLY_FIRST_NOCI_ELEMENTS) then call NonOrthogonalCI_diagonalizeCImatrix(NonOrthogonalCI_instance) - call NonOrthogonalCI_generateDensities(NonOrthogonalCI_instance) + call NonOrthogonalCI_generateSuperposedSystem(NonOrthogonalCI_instance) + call NonOrthogonalCI_buildDensityMatrix(NonOrthogonalCI_instance) + call NonOrthogonalCI_getNaturalOrbitals(NonOrthogonalCI_instance) + call NonOrthogonalCI_computeFranckCondon(NonOrthogonalCI_instance) + call NonOrthogonalCI_saveToFile(NonOrthogonalCI_instance) else write (*,"(T10,A)") "COMPUTED NOCI ELEMENTS ONLY WITH RESPECT TO THE FIRST GEOMETRY - YOU HAVE TO SOLVE THE CI EQUATION MANUALLY!" end if @@ -206,7 +210,11 @@ program NOCI call NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(NonOrthogonalCI_instance) if(.not. CONTROL_instance%ONLY_FIRST_NOCI_ELEMENTS) then call NonOrthogonalCI_diagonalizeCImatrix(NonOrthogonalCI_instance) - call NonOrthogonalCI_generateDensities(NonOrthogonalCI_instance) + call NonOrthogonalCI_generateSuperposedSystem(NonOrthogonalCI_instance) + call NonOrthogonalCI_buildDensityMatrix(NonOrthogonalCI_instance) + call NonOrthogonalCI_getNaturalOrbitals(NonOrthogonalCI_instance) + call NonOrthogonalCI_computeFranckCondon(NonOrthogonalCI_instance) + call NonOrthogonalCI_saveToFile(NonOrthogonalCI_instance) else write (*,"(T10,A)") "COMPUTED NOCI ELEMENTS ONLY WITH RESPECT TO THE FIRST GEOMETRY - YOU HAVE TO SOLVE THE CI EQUATION MANUALLY!" end if diff --git a/src/NOCI/NonOrthogonalCI.f90 b/src/NOCI/NonOrthogonalCI.f90 index b9642341..756d2f01 100644 --- a/src/NOCI/NonOrthogonalCI.f90 +++ b/src/NOCI/NonOrthogonalCI.f90 @@ -70,24 +70,27 @@ module NonOrthogonalCI_ integer :: numberOfEquivalentSystems integer :: numberOfTransformedCenters integer :: numberOfIndividualTransformations - integer :: numberOfUniqueSystems !sort of symmetry - integer :: numberOfUniquePairs !sort of symmetry integer :: printMatrixThreshold - type(IVector) :: systemTypes !sort of symmetry integer, allocatable :: rotationCenterList(:,:) type(Matrix) :: configurationOverlapMatrix, configurationHamiltonianMatrix, configurationCoefficients - type(Matrix), allocatable :: configurationKineticMatrix(:), configurationPuntualMatrix(:), configurationExternalMatrix(:) - type(Matrix), allocatable :: configurationTwoParticlesMatrix(:,:) - type(IMatrix) :: configurationPairTypes !, uniqueOverlapElements, uniqueHamiltonianElements + type(Matrix), allocatable :: configurationKineticMatrix(:), configurationPuntualMatrix(:), configurationExternalMatrix(:), configurationExchangeMatrix(:) + type(Matrix), allocatable :: configurationHartreeMatrix(:,:) type(Vector) :: configurationCorrelationEnergies, statesEigenvalues + type(IVector), allocatable :: sysBasisList(:,:) type(Matrix), allocatable :: HFCoefficients(:,:) - type(MolecularSystem), allocatable :: MolecularSystems(:) - type(MolecularSystem), allocatable :: uniqueMolecularSystems(:) - ! type(Matrix), allocatable :: HCoreMatrices(:,:,:) - ! type(Matrix), allocatable :: inverseOverlapMatrices(:,:,:) + type(Matrix), allocatable :: mergedCoefficients(:) + type(Matrix), allocatable :: mergedOverlapMatrix(:) + type(Matrix), allocatable :: mergedDensityMatrix(:,:) + type(MolecularSystem), allocatable :: molecularSystems(:) + type(MolecularSystem) :: mergedMolecularSystem character(50) :: transformationType character(15),allocatable :: systemLabels(:) real(8) :: refEnergy + ! integer :: numberOfUniqueSystems !sort of symmetry + ! integer :: numberOfUniquePairs !sort of symmetry + ! type(IVector) :: systemTypes !sort of symmetry + ! type(IMatrix) :: configurationPairTypes !, uniqueOverlapElements, uniqueHamiltonianElements + ! type(MolecularSystem), allocatable :: uniqueMolecularSystems(:) end type NonOrthogonalCI type(NonOrthogonalCI), public :: NonOrthogonalCI_instance @@ -99,7 +102,11 @@ module NonOrthogonalCI_ NonOrthogonalCI_runHFs,& NonOrthogonalCI_buildOverlapAndHamiltonianMatrix,& NonOrthogonalCI_diagonalizeCImatrix,& - NonOrthogonalCI_generateDensities + NonOrthogonalCI_generateSuperposedSystem,& + NonOrthogonalCI_buildDensityMatrix,& + NonOrthogonalCI_getNaturalOrbitals,& + NonOrthogonalCI_saveToFile,& + NonOrthogonalCI_computeFranckCondon private @@ -126,9 +133,9 @@ subroutine NonOrthogonalCI_constructor(this) this%numberOfEllipsoidRejectedSystems=0 this%numberOfPPdistanceRejectedSystems=0 this%numberOfNPdistanceRejectedSystems=0 - this%numberOfUniqueSystems=0 - this%numberOfUniquePairs=0 - this%printMatrixThreshold=10 + ! this%numberOfUniqueSystems=0 + ! this%numberOfUniquePairs=0 + this%printMatrixThreshold=30 numberOfTranslationCenters=0 numberOfRotationCenters=0 @@ -215,10 +222,14 @@ subroutine NonOrthogonalCI_constructor(this) ! call Vector_constructorInteger(this%systemTypes,this%numberOfIndividualTransformations**this%numberOfTransformedCenters,0) - allocate(this%configurationKineticMatrix(molecularSystem_instance%numberOfQuantumSpecies),& + allocate(this%mergedDensityMatrix(CONTROL_instance%CI_STATES_TO_PRINT,molecularSystem_instance%numberOfQuantumSpecies),& + this%mergedOverlapMatrix(molecularSystem_instance%numberOfQuantumSpecies),& + this%mergedCoefficients(molecularSystem_instance%numberOfQuantumSpecies),& + this%configurationKineticMatrix(molecularSystem_instance%numberOfQuantumSpecies),& this%configurationPuntualMatrix(molecularSystem_instance%numberOfQuantumSpecies),& this%configurationExternalMatrix(molecularSystem_instance%numberOfQuantumSpecies),& - this%configurationTwoParticlesMatrix(molecularSystem_instance%numberOfQuantumSpecies,molecularSystem_instance%numberOfQuantumSpecies)) + this%configurationExchangeMatrix(molecularSystem_instance%numberOfQuantumSpecies),& + this%configurationHartreeMatrix(molecularSystem_instance%numberOfQuantumSpecies,molecularSystem_instance%numberOfQuantumSpecies)) end subroutine NonOrthogonalCI_constructor !> @@ -238,8 +249,6 @@ subroutine NonOrthogonalCI_displaceGeometries(this) integer :: coordsUnit integer :: i,j integer :: closestSystem - integer :: systemType - logical :: newSystemFlag logical :: skip real(8) :: timeA @@ -248,7 +257,7 @@ subroutine NonOrthogonalCI_displaceGeometries(this) call MolecularSystem_copyConstructor(originalMolecularSystem, molecularSystem_instance) !!Dynamically allocated through the displacement routine - allocate(this%MolecularSystems(0)) + allocate(this%molecularSystems(0)) allocate(transformationCounter(this%numberOfTransformedCenters)) @@ -292,8 +301,8 @@ subroutine NonOrthogonalCI_displaceGeometries(this) call MolecularSystem_showCartesianMatrix(displacedMolecularSystem,unit=coordsUnit) !Classify the system according to its distance matrix (symmetry) - if(CONTROL_instance%CONFIGURATION_USE_SYMMETRY .eqv. .true.) & - call NonOrthogonalCI_classifyNewSystem(this,systemType, newSystemFlag) + ! if(CONTROL_instance%CONFIGURATION_USE_SYMMETRY .eqv. .true.) & + ! call NonOrthogonalCI_classifyNewSystem(this,systemType, newSystemFlag) !Check if the new system is not beyond the max displacement if(skip) then @@ -372,10 +381,8 @@ subroutine NonOrthogonalCI_displaceGeometries(this) print *, "" - if(CONTROL_instance%CONFIGURATION_USE_SYMMETRY .eqv. .true.) & - call Matrix_constructorInteger(this%configurationPairTypes,int(this%numberOfDisplacedSystems,8), int(this%numberOfDisplacedSystems,8), 0) - - + ! if(CONTROL_instance%CONFIGURATION_USE_SYMMETRY .eqv. .true.) & + ! call Matrix_constructorInteger(this%configurationPairTypes,int(this%numberOfDisplacedSystems,8), int(this%numberOfDisplacedSystems,8), 0) ! minEnergy=0.0 !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time displacing coordinates : ", omp_get_wtime() - timeA ," (s)" @@ -419,7 +426,7 @@ subroutine NonOrthogonalCI_readGeometries(this) read(coordsUnit,*) string, this%numberOfDisplacedSystems print *, "reading ", this%numberOfDisplacedSystems, " systems" - allocate(this%MolecularSystems(this%numberOfDisplacedSystems)) + allocate(this%molecularSystems(this%numberOfDisplacedSystems)) do sysI = 1, this%numberOfDisplacedSystems call MolecularSystem_copyConstructor(molecularSystem_instance, originalMolecularSystem) @@ -463,7 +470,7 @@ subroutine NonOrthogonalCI_readGeometries(this) molecularSystem_instance%pointCharges(i)%origin = origin/AMSTRONG end do - call MolecularSystem_copyConstructor(this%MolecularSystems(sysI), molecularSystem_instance) + call MolecularSystem_copyConstructor(this%molecularSystems(sysI), molecularSystem_instance) end do @@ -653,7 +660,7 @@ subroutine NonOrthogonalCI_checkNewSystemDisplacement(this,newMolecularSystem,cl do sysI=1, this%numberOfDisplacedSystems - call MolecularSystem_GetTwoSystemsDisplacement(this%MolecularSystems(sysI), newMolecularSystem, displacementVector) + call MolecularSystem_GetTwoSystemsDisplacement(this%molecularSystems(sysI), newMolecularSystem, displacementVector) dispSum=0.0 do i=1, newMolecularSystem%numberOfQuantumSpecies @@ -751,65 +758,65 @@ end subroutine NonOrthogonalCI_checkSameChargesDistance !! !! @param this, systemType: integer defining system equivalence type, newSystemFlag: returns if the system is new or not !< - subroutine NonOrthogonalCI_classifyNewSystem(this, systemType, newSystemFlag) - implicit none - type(NonOrthogonalCI) :: this - integer :: systemType - logical :: newSystemFlag + ! subroutine NonOrthogonalCI_classifyNewSystem(this, systemType, newSystemFlag) + ! implicit none + ! type(NonOrthogonalCI) :: this + ! integer :: systemType + ! logical :: newSystemFlag - type(MolecularSystem) :: currentMolecularSystem - type(Matrix) :: currentDistanceMatrix,previousDistanceMatrix + ! type(MolecularSystem) :: currentMolecularSystem + ! type(Matrix) :: currentDistanceMatrix,previousDistanceMatrix - integer :: sysI, i, checkingType - logical :: match + ! integer :: sysI, i, checkingType + ! logical :: match - call MolecularSystem_copyConstructor(currentMolecularSystem, molecularSystem_instance) - systemType=0 - newSystemFlag=.true. - currentDistanceMatrix=ParticleManager_getDistanceMatrix() + ! call MolecularSystem_copyConstructor(currentMolecularSystem, molecularSystem_instance) + ! systemType=0 + ! newSystemFlag=.true. + ! currentDistanceMatrix=ParticleManager_getDistanceMatrix() - ! print *, "Current distance matrix" - ! call Matrix_show(currentDistanceMatrix) + ! ! print *, "Current distance matrix" + ! ! call Matrix_show(currentDistanceMatrix) - types: do checkingType=1, this%numberOfUniqueSystems - ! print *, "checkingType", checkingType - systems: do sysI=1, this%numberOfDisplacedSystems + ! types: do checkingType=1, this%numberOfUniqueSystems + ! ! print *, "checkingType", checkingType + ! systems: do sysI=1, this%numberOfDisplacedSystems - if(this%systemTypes%values(sysI) .eq. checkingType) then - call MolecularSystem_copyConstructor(molecularSystem_instance, this%MolecularSystems(sysI)) + ! if(this%systemTypes%values(sysI) .eq. checkingType) then + ! call MolecularSystem_copyConstructor(molecularSystem_instance, this%molecularSystems(sysI)) - previousDistanceMatrix=ParticleManager_getDistanceMatrix() + ! previousDistanceMatrix=ParticleManager_getDistanceMatrix() - ! print *, "Comparing with previous distance matrix", checkingType - ! call Matrix_show(previousDistanceMatrix) + ! ! print *, "Comparing with previous distance matrix", checkingType + ! ! call Matrix_show(previousDistanceMatrix) - match=.true. - do i=1, size(currentDistanceMatrix%values(:,1)) - if(sum(abs(currentDistanceMatrix%values(i,:) - previousDistanceMatrix%values(i,:))) .gt. & - CONTROL_instance%CONFIGURATION_EQUIVALENCE_DISTANCE ) then - match=.false. - exit - end if - end do + ! match=.true. + ! do i=1, size(currentDistanceMatrix%values(:,1)) + ! if(sum(abs(currentDistanceMatrix%values(i,:) - previousDistanceMatrix%values(i,:))) .gt. & + ! CONTROL_instance%CONFIGURATION_EQUIVALENCE_DISTANCE ) then + ! match=.false. + ! exit + ! end if + ! end do - ! print *, "match?", match - - if(match) then - systemType=this%systemTypes%values(sysI) - newSystemFlag=.false. - exit types - else - cycle types - end if - end if - end do systems - end do types + ! ! print *, "match?", match + + ! if(match) then + ! systemType=this%systemTypes%values(sysI) + ! newSystemFlag=.false. + ! exit types + ! else + ! cycle types + ! end if + ! end if + ! end do systems + ! end do types - ! print *, "newSystemFlag", newSystemFlag + ! ! print *, "newSystemFlag", newSystemFlag - call MolecularSystem_copyConstructor(molecularSystem_instance, currentMolecularSystem) + ! call MolecularSystem_copyConstructor(molecularSystem_instance, currentMolecularSystem) - end subroutine NonOrthogonalCI_classifyNewSystem + ! end subroutine NonOrthogonalCI_classifyNewSystem !> !! @brief Run a Hartree-Fock calculation at displaced geometries and fill CI matrix diagonals @@ -842,8 +849,10 @@ subroutine NonOrthogonalCI_runHFs(this) int(this%numberOfDisplacedSystems,8), int(this%numberOfDisplacedSystems,8), 0.0_8) call Matrix_constructor(this%configurationExternalMatrix(speciesID), & int(this%numberOfDisplacedSystems,8), int(this%numberOfDisplacedSystems,8), 0.0_8) + call Matrix_constructor(this%configurationExchangeMatrix(speciesID), & + int(this%numberOfDisplacedSystems,8), int(this%numberOfDisplacedSystems,8), 0.0_8) do otherSpeciesID=speciesID, MolecularSystem_instance%numberOfQuantumSpecies - call Matrix_constructor(this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID), & + call Matrix_constructor(this%configurationHartreeMatrix(speciesID,otherSpeciesID), & int(this%numberOfDisplacedSystems,8), int(this%numberOfDisplacedSystems,8), 0.0_8) end do end do @@ -862,10 +871,10 @@ subroutine NonOrthogonalCI_runHFs(this) write (coordsUnit,'(A25,I20)') "numberOfDisplacedSystems ", this%numberOfDisplacedSystems do sysI=1, this%numberOfDisplacedSystems - write(this%systemLabels(sysI), '(A)') trim(this%MolecularSystems(sysI)%description) + write(this%systemLabels(sysI), '(A)') trim(this%molecularSystems(sysI)%description) !!Do SCF without calling lowdin-scf.x - call MolecularSystem_copyConstructor(molecularSystem_instance, this%MolecularSystems(sysI)) + call MolecularSystem_copyConstructor(molecularSystem_instance, this%molecularSystems(sysI)) CONTROL_instance%PRINT_LEVEL=0 if ( CONTROL_instance%METHOD .eq. "RKS" .or. CONTROL_instance%METHOD .eq. "UKS" ) & @@ -906,6 +915,8 @@ subroutine NonOrthogonalCI_runHFs(this) call MultiSCF_solveHartreeFockRoothan(MultiSCF_instance,WaveFunction_instance,Libint2Instance) !Save HF results + ! call MultiSCF_saveWfn(MultiSCF_instance,WaveFunction_instance) + call MolecularSystem_copyConstructor(this%molecularSystems(sysI),molecularSystem_instance) this%configurationHamiltonianMatrix%values(sysI,sysI)=MultiSCF_instance%totalEnergy do speciesID = 1, molecularSystem_instance%numberOfQuantumSpecies @@ -914,11 +925,9 @@ subroutine NonOrthogonalCI_runHFs(this) this%configurationKineticMatrix(speciesID)%values(sysI,sysI)=WaveFunction_instance(speciesID)%kineticEnergy this%configurationPuntualMatrix(speciesID)%values(sysI,sysI)=WaveFunction_instance(speciesID)%puntualInteractionEnergy this%configurationExternalMatrix(speciesID)%values(sysI,sysI)=WaveFunction_instance(speciesID)%externalPotentialEnergy - this%configurationTwoParticlesMatrix(speciesID,speciesID)%values(sysI,sysI)=& - WaveFunction_instance(speciesID)%hartreeEnergy(speciesID)+WaveFunction_instance(speciesID)%exchangeHFEnergy - - do otherSpeciesID = speciesID+1, molecularSystem_instance%numberOfQuantumSpecies - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(sysI,sysI)=& + this%configurationExchangeMatrix(speciesID)%values(sysI,sysI)=WaveFunction_instance(speciesID)%exchangeHFEnergy + do otherSpeciesID = speciesID, molecularSystem_instance%numberOfQuantumSpecies + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,sysI)=& WaveFunction_instance(speciesID)%hartreeEnergy(otherSpeciesID) end do end do @@ -955,6 +964,9 @@ subroutine NonOrthogonalCI_runHFs(this) write (*,'(A10,I10,A10,ES20.12,A20,ES20.12)') "Geometry ", sysI, "Energy", this%configurationHamiltonianMatrix%values(sysI,sysI), & "Correlation energy", this%configurationCorrelationEnergies%values(sysI) call MolecularSystem_showCartesianMatrix(MolecularSystem_instance) + ! do speciesID = 1, MolecularSystem_instance%numberOfQuantumSpecies + ! print *, "sysI", sysI, "speciesID", speciesID, "occupation number", MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) + ! end do end if call DirectIntegralManager_destructor(Libint2Instance) @@ -1023,7 +1035,7 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) type(NonOrthogonalCI) :: this type(MolecularSystem), allocatable :: mergedMolecularSystem(:) type(Libint2Interface), allocatable :: Libint2ParallelInstance(:,:) - integer, allocatable :: sysIIbatch(:) + integer, allocatable :: sysIbatch(:), sysIIbatch(:) integer :: sysI,sysII,me,mySysII type(Matrix), allocatable :: mergedCoefficients(:), inverseOverlapMatrices(:) type(IVector), allocatable :: sysIbasisList(:,:),sysIIbasisList(:,:) @@ -1071,14 +1083,12 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) ncores=CONTROL_instance%NUMBER_OF_CORES batchSize=this%numberOfDisplacedSystems print *, "ncores", ncores, "batchsize", batchSize - if (this%numberOfDisplacedSystems .le. this%printMatrixThreshold) then - write (*,'(A10,A10,A20,A20)') "sysI", "sysII", "OverlapElement", "HamiltonianElement" - end if allocate(mergedMolecularSystem(batchSize),& mergedCoefficients(nspecies),& inverseOverlapMatrices(nspecies),& Libint2ParallelInstance(nspecies,batchSize),& + sysIbatch(batchSize),& sysIIbatch(batchSize),& sysIbasisList(nspecies,batchSize),& sysIIbasisList(nspecies,batchSize)) @@ -1128,7 +1138,7 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) !This generates a new molecular system ! print *, "Merging systems from geometries ", sysI, mySysII call MolecularSystem_mergeTwoSystems(mergedMolecularSystem(me), & - this%MolecularSystems(sysI), this%MolecularSystems(mySysII),sysIbasisList(1:nspecies,me),sysIIbasisList(1:nspecies,me)) + this%molecularSystems(sysI), this%molecularSystems(mySysII),sysIbasisList(1:nspecies,me),sysIIbasisList(1:nspecies,me)) ! call MolecularSystem_showInformation() ! call MolecularSystem_showParticlesInformation() ! call MolecularSystem_showCartesianMatrix(mergedMolecularSystem) @@ -1155,7 +1165,7 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) ! cycle systemII !! Merge occupied coefficients into a single matrix call NonOrthogonalCI_mergeCoefficients(this%HFCoefficients(sysI,:),this%HFCoefficients(mySysII,:),& - this%MolecularSystems(sysI),this%MolecularSystems(mySysII),mergedMolecularSystem(me),& + this%molecularSystems(sysI),this%molecularSystems(mySysII),mergedMolecularSystem(me),& sysIbasisList(1:nspecies,me),sysIIbasisList(1:nspecies,me),mergedCoefficients) !$ timeA = omp_get_wtime() @@ -1202,19 +1212,27 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) end if !Yu2020 magical empirical correction - if(CONTROL_instance%EMPIRICAL_OVERLAP_CORRECTION .and. abs(this%configurationOverlapMatrix%values(sysI,mySysII)) .gt. CONTROL_instance%CONFIGURATION_OVERLAP_THRESHOLD) then + if(CONTROL_instance%EMPIRICAL_OVERLAP_CORRECTION .and. & + abs(this%configurationOverlapMatrix%values(sysI,mySysII)) .gt. CONTROL_instance%CONFIGURATION_OVERLAP_THRESHOLD) then empiricalScaleFactor=CONTROL_instance%EMPIRICAL_OVERLAP_PARAMETER_A*& abs(this%configurationOverlapMatrix%values(sysI,mySysII))**CONTROL_instance%EMPIRICAL_OVERLAP_PARAMETER_B/& abs(this%configurationOverlapMatrix%values(sysI,mySysII)) - this%configurationOverlapMatrix%values(sysI,mySysII)=this%configurationOverlapMatrix%values(sysI,mySysII)*empiricalScaleFactor - this%configurationHamiltonianMatrix%values(sysI,mySysII)=this%configurationHamiltonianMatrix%values(sysI,mySysII)*empiricalScaleFactor + this%configurationOverlapMatrix%values(sysI,mySysII)=& + this%configurationOverlapMatrix%values(sysI,mySysII)*empiricalScaleFactor + this%configurationHamiltonianMatrix%values(sysI,mySysII)=& + this%configurationHamiltonianMatrix%values(sysI,mySysII)*empiricalScaleFactor do speciesID=1, nspecies - this%configurationKineticMatrix(speciesID)%values(sysI,mySysII)=this%configurationKineticMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor - this%configurationPuntualMatrix(speciesID)%values(sysI,mySysII)=this%configurationPuntualMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor - this%configurationExternalMatrix(speciesID)%values(sysI,mySysII)=this%configurationExternalMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor + this%configurationKineticMatrix(speciesID)%values(sysI,mySysII)=& + this%configurationKineticMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor + this%configurationPuntualMatrix(speciesID)%values(sysI,mySysII)=& + this%configurationPuntualMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor + this%configurationExternalMatrix(speciesID)%values(sysI,mySysII)=& + this%configurationExternalMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor + this%configurationExchangeMatrix(speciesID)%values(sysI,mySysII)=& + this%configurationExchangeMatrix(speciesID)%values(sysI,mySysII)*empiricalScaleFactor do otherSpeciesID=speciesID, nspecies - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII)=& - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII) + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII)=& + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII)*empiricalScaleFactor end do end do end if @@ -1227,17 +1245,40 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) this%configurationKineticMatrix(speciesID)%values(mySysII,sysI)=this%configurationKineticMatrix(speciesID)%values(sysI,mySysII) this%configurationPuntualMatrix(speciesID)%values(mySysII,sysI)=this%configurationPuntualMatrix(speciesID)%values(sysI,mySysII) this%configurationExternalMatrix(speciesID)%values(mySysII,sysI)=this%configurationExternalMatrix(speciesID)%values(sysI,mySysII) + this%configurationExchangeMatrix(speciesID)%values(mySysII,sysI)=this%configurationExchangeMatrix(speciesID)%values(sysI,mySysII) do otherSpeciesID=speciesID, nspecies - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(mySysII,sysI)=& - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII) + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(mySysII,sysI)=& + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII) end do end do write (matrixUnit,'(I10,I10,ES20.12,ES20.12)') sysI, mySysII, & this%configurationOverlapMatrix%values(sysI,mySysII), this%configurationHamiltonianMatrix%values(sysI,mySysII) + if (this%numberOfDisplacedSystems .le. this%printMatrixThreshold) then - write (*,'(I10,I10,ES20.12,ES20.12)') sysI, mySysII, & - this%configurationOverlapMatrix%values(sysI,mySysII), this%configurationHamiltonianMatrix%values(sysI,mySysII) + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, "Overlap element = ", this%configurationOverlapMatrix%values(sysI,mySysII) + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, "Hamiltonian element = ", this%configurationHamiltonianMatrix%values(sysI,mySysII) + + do speciesID = 1, nspecies + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, trim( this%molecularSystems(sysI)%species(speciesID)%name ) // & + " Kinetic element = ", this%configurationKineticMatrix(speciesID)%values(sysI,mySysII) + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, trim( this%molecularSystems(sysI)%species(speciesID)%name ) // & + " Puntual element = ", this%configurationPuntualMatrix(speciesID)%values(sysI,mySysII) + if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, trim( this%molecularSystems(sysI)%species(speciesID)%name ) // & + " External element = ", this%configurationExternalMatrix(speciesID)%values(sysI,mySysII) + end do + do speciesID=1, nspecies + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, trim( this%molecularSystems(sysI)%species(speciesID)%name ) // & + " Exchange element = ", this%configurationExchangeMatrix(speciesID)%values(sysI,mySysII) + do otherSpeciesID=speciesID, nspecies + write (*,'(I10,I10,A38,ES20.12)') sysI, mySysII, trim( MolecularSystem_instance%species(speciesID)%name ) // & + "/"//trim( MolecularSystem_instance%species(otherSpeciesID)%name ) // & + " Hartree element = ", this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,mySysII) + end do + end do + print *, "" + end if call DirectIntegralManager_destructor(Libint2ParallelInstance(1:nspecies,me)) @@ -1270,6 +1311,7 @@ subroutine NonOrthogonalCI_buildOverlapAndHamiltonianMatrix(this) mergedCoefficients,& inverseOverlapMatrices,& Libint2ParallelInstance,& + sysIbatch,& sysIIbatch,& sysIbasisList,& sysIIbasisList) @@ -1452,25 +1494,25 @@ subroutine NonOrthogonalCI_prescreenOverlap(this,sysI,sysII,estimatedOverlap) real(8) :: massThreshold, minExponent, speciesOverlap !displacement vectors contains the max distance between equivalent basis function centers - allocate(displacementVector(this%MolecularSystems(sysI)%numberOfQuantumSpecies)) + allocate(displacementVector(this%molecularSystems(sysI)%numberOfQuantumSpecies)) - call MolecularSystem_GetTwoSystemsDisplacement(this%MolecularSystems(sysI), this%MolecularSystems(sysII),displacementVector(:)) + call MolecularSystem_GetTwoSystemsDisplacement(this%molecularSystems(sysI), this%molecularSystems(sysII),displacementVector(:)) estimatedOverlap=1.0 !only compute for heavy particles, maybe should be a control parameter massThreshold=10.0 - do speciesID = 1, this%MolecularSystems(sysI)%numberOfQuantumSpecies - if(this%MolecularSystems(sysI)%species(speciesID)%mass .lt. massThreshold) cycle + do speciesID = 1, this%molecularSystems(sysI)%numberOfQuantumSpecies + if(this%molecularSystems(sysI)%species(speciesID)%mass .lt. massThreshold) cycle speciesOverlap=1.0 !!get smallest exponent of the basis set - do k = 1, size(this%MolecularSystems(sysI)%species(speciesID)%particles) + do k = 1, size(this%molecularSystems(sysI)%species(speciesID)%particles) minExponent=1.0E8 - do l = 1, size(this%MolecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction) - do m = 1, size(this%MolecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction(l)%orbitalExponents) - if(this%MolecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction(l)%orbitalExponents(m).lt.minExponent) & - minExponent=this%MolecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction(l)%orbitalExponents(m) + do l = 1, size(this%molecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction) + do m = 1, size(this%molecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction(l)%orbitalExponents) + if(this%molecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction(l)%orbitalExponents(m).lt.minExponent) & + minExponent=this%molecularSystems(sysI)%species(speciesID)%particles(k)%basis%contraction(l)%orbitalExponents(m) !Assume a 1S GTF ! normCoefficients(speciesID)=(2.0*minExponents(speciesID)/Math_PI)**(3.0/4.0) end do @@ -1605,7 +1647,7 @@ subroutine NonOrthogonalCI_computeOverlapAndHCoreElements(this,sysI,sysII,merged do speciesID = 1, MolecularSystem_instance%numberOfQuantumSpecies numberOfContractions=MolecularSystem_getTotalNumberOfContractions(speciesID,mergedMolecularSystem) - occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) particlesPerOrbital=MolecularSystem_getEta(speciesID,mergedMolecularSystem) !! Calculate one- particle integrals call DirectIntegralManager_getOverlapIntegrals(mergedMolecularSystem,speciesID,& @@ -1671,7 +1713,7 @@ subroutine NonOrthogonalCI_computeOverlapAndHCoreElements(this,sysI,sysII,merged do speciesID = 1, MolecularSystem_instance%numberOfQuantumSpecies numberOfContractions=MolecularSystem_getTotalNumberOfContractions(speciesID,mergedMolecularSystem) - occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) particlesPerOrbital=MolecularSystem_getEta(speciesID,mergedMolecularSystem) call Matrix_constructor(auxKineticMatrix(speciesID),& @@ -1749,7 +1791,7 @@ subroutine NonOrthogonalCI_computeOverlapAndHCoreElements(this,sysI,sysII,merged oneParticleKineticEnergy=0.0 oneParticleAttractionEnergy=0.0 oneParticleExternalEnergy=0.0 - occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) do a=1, occupationNumber !sysI do b=1, occupationNumber !sysII oneParticleKineticEnergy=oneParticleKineticEnergy+ molecularKineticMatrix(speciesID)%values(a,b)*& @@ -1793,7 +1835,7 @@ subroutine NonOrthogonalCI_twoParticlesContributions(this,sysI,sysII,mergedMolec integer :: otherNumberOfContractions,otherOccupationNumber,otherParticlesPerOrbital integer :: ssize1, auxIndex, auxIndex1 integer :: a,b,bb,c,d,dd,i,j - real(8) :: interactionEnergy + real(8) :: hartreeEnergy, exchangeEnergy allocate(fourCenterIntegrals(mergedMolecularSystem%numberOfQuantumSpecies,mergedMolecularSystem%numberOfQuantumSpecies), & twoIndexArray(mergedMolecularSystem%numberOfQuantumSpecies), & @@ -1858,9 +1900,10 @@ subroutine NonOrthogonalCI_twoParticlesContributions(this,sysI,sysII,mergedMolec !!Same species repulsion do i=1, mergedMolecularSystem%numberOfQuantumSpecies numberOfContractions=MolecularSystem_getTotalNumberOfContractions(i,mergedMolecularSystem) - occupationNumber=MolecularSystem_getOcupationNumber(i,this%MolecularSystems(sysI)) + occupationNumber=MolecularSystem_getOcupationNumber(i,this%molecularSystems(sysI)) particlesPerOrbital=MolecularSystem_getEta(i,mergedMolecularSystem) - interactionEnergy=0.0 + hartreeEnergy=0.0 + exchangeEnergy=0.0 do a=1,occupationNumber !sysI do b=occupationNumber+1, 2*occupationNumber !sysII bb=b-occupationNumber @@ -1868,43 +1911,45 @@ subroutine NonOrthogonalCI_twoParticlesContributions(this,sysI,sysII,mergedMolec do d=occupationNumber+1, 2*occupationNumber !sysII dd=d-occupationNumber auxIndex = fourIndexArray(i)%values(twoIndexArray(i)%values(a,b), twoIndexArray(i)%values(c,d) ) - interactionEnergy=interactionEnergy+0.5*fourCenterIntegrals(i,i)%values(auxIndex, 1)*& - (inverseOverlapMatrices(i)%values(bb,a)*inverseOverlapMatrices(i)%values(dd,c)*particlesPerOrbital**2& !coulomb - -inverseOverlapMatrices(i)%values(dd,a)*inverseOverlapMatrices(i)%values(bb,c)*particlesPerOrbital) !exchange + hartreeEnergy=hartreeEnergy+0.5*fourCenterIntegrals(i,i)%values(auxIndex, 1)*& + inverseOverlapMatrices(i)%values(bb,a)*inverseOverlapMatrices(i)%values(dd,c)*particlesPerOrbital**2 !coulomb + exchangeEnergy=exchangeEnergy-0.5*fourCenterIntegrals(i,i)%values(auxIndex, 1)*& + inverseOverlapMatrices(i)%values(dd,a)*inverseOverlapMatrices(i)%values(bb,c)*particlesPerOrbital !exchange ! print *, a, b, c, d, twoIndexArray(i)%values(a,b), twoIndexArray(i)%values(c,d), fourIndexArray(i)%values( & ! twoIndexArray(i)%values(a,b), twoIndexArray(i)%values(c,d)), end do end do end do end do - this%configurationTwoParticlesMatrix(i,i)%values(sysI,sysII)=interactionEnergy*this%configurationOverlapMatrix%values(sysI,sysII) + this%configurationHartreeMatrix(i,i)%values(sysI,sysII)=hartreeEnergy*this%configurationOverlapMatrix%values(sysI,sysII) + this%configurationExchangeMatrix(i)%values(sysI,sysII)=exchangeEnergy*this%configurationOverlapMatrix%values(sysI,sysII) this%configurationHamiltonianMatrix%values(sysI,sysII)=this%configurationHamiltonianMatrix%values(sysI,sysII)+& - interactionEnergy*this%configurationOverlapMatrix%values(sysI,sysII) - ! print *, "same species interactionEnergy for species", i, interactionEnergy + (hartreeEnergy+exchangeEnergy)*this%configurationOverlapMatrix%values(sysI,sysII) + ! print *, "same species interactionEnergy for species", i, hartreeEnergy, exchangeEnergy end do !!Interspecies repulsion do i=1, mergedMolecularSystem%numberOfQuantumSpecies-1 numberOfContractions=MolecularSystem_getTotalNumberOfContractions(i,mergedMolecularSystem) - occupationNumber=MolecularSystem_getOcupationNumber(i,this%MolecularSystems(sysI)) + occupationNumber=MolecularSystem_getOcupationNumber(i,this%molecularSystems(sysI)) particlesPerOrbital=MolecularSystem_getEta(i,mergedMolecularSystem) do j=i+1, mergedMolecularSystem%numberOfQuantumSpecies otherNumberOfContractions=MolecularSystem_getTotalNumberOfContractions(j,mergedMolecularSystem) otherOccupationNumber=MolecularSystem_getOcupationNumber(j,mergedMolecularSystem) otherParticlesPerOrbital=MolecularSystem_getEta(j,mergedMolecularSystem) - interactionEnergy=0.0 + hartreeEnergy=0.0 ssize1 = max(otherNumberOfContractions,otherOccupationNumber) ssize1 = ( ssize1 * ( ssize1 + 1 ) ) / 2 - otherOccupationNumber=MolecularSystem_getOcupationNumber(j,this%MolecularSystems(sysI)) + otherOccupationNumber=MolecularSystem_getOcupationNumber(j,this%molecularSystems(sysI)) do a=1, occupationNumber !sysI do b=occupationNumber+1, 2*occupationNumber !sysII - bb=b-MolecularSystem_getOcupationNumber(i,this%MolecularSystems(sysI)) + bb=b-MolecularSystem_getOcupationNumber(i,this%molecularSystems(sysI)) auxIndex1 = ssize1 * (twoIndexArray(i)%values(a,b) - 1 ) do c=1, otherOccupationNumber !sysI do d=otherOccupationNumber+1,2*otherOccupationNumber !sysII dd=d-otherOccupationNumber auxIndex = auxIndex1 + twoIndexArray(j)%values(c,d) - interactionEnergy=interactionEnergy+fourCenterIntegrals(i,j)%values(auxIndex, 1)*& + hartreeEnergy=hartreeEnergy+fourCenterIntegrals(i,j)%values(auxIndex, 1)*& inverseOverlapMatrices(i)%values(bb,a)*inverseOverlapMatrices(j)%values(dd,c)*& particlesPerOrbital*otherParticlesPerOrbital ! print *, a, b, c, d, fourCenterIntegrals(i,j)%values(auxIndex, 1), inverseOverlapMatrices(i)%values(bb,a), inverseOverlapMatrices(j)%values(dd,c) @@ -1912,10 +1957,10 @@ subroutine NonOrthogonalCI_twoParticlesContributions(this,sysI,sysII,mergedMolec end do end do end do - this%configurationTwoParticlesMatrix(i,j)%values(sysI,sysII)=interactionEnergy*this%configurationOverlapMatrix%values(sysI,sysII) + this%configurationHartreeMatrix(i,j)%values(sysI,sysII)=hartreeEnergy*this%configurationOverlapMatrix%values(sysI,sysII) this%configurationHamiltonianMatrix%values(sysI,sysII)=this%configurationHamiltonianMatrix%values(sysI,sysII)+& - interactionEnergy*this%configurationOverlapMatrix%values(sysI,sysII) - ! print *, "interspecies interactionEnergy for species", i, j, interactionEnergy + hartreeEnergy*this%configurationOverlapMatrix%values(sysI,sysII) + ! print *, "interspecies hartreeEnergy for species", i, j, hartreeEnergy end do end do @@ -1923,6 +1968,11 @@ subroutine NonOrthogonalCI_twoParticlesContributions(this,sysI,sysII,mergedMolec end subroutine NonOrthogonalCI_twoParticlesContributions + !> + !! @brief Solves the NOCI matrix equation + !! + !! @param this + !< subroutine NonOrthogonalCI_diagonalizeCImatrix(this) implicit none type(NonOrthogonalCI) :: this @@ -2023,7 +2073,7 @@ subroutine NonOrthogonalCI_diagonalizeCImatrix(this) write(*,"(A)") " EIGENVALUES AND EIGENVECTORS: " write(*,"(A)") "=========================================" write(*,"(A)") "" - do state = 1, CONTROL_instance%NUMBER_OF_CI_STATES + do state = 1, min(CONTROL_instance%NUMBER_OF_CI_STATES,this%numberOfDisplacedSystems) write (*,"(A)") "" write (*,"(T9,A17,I3,A10, F25.12)") "STATE: ", state, " ENERGY = ", this%statesEigenvalues%values(state) write (*,"(A38)") "Components: " @@ -2079,22 +2129,37 @@ subroutine NonOrthogonalCI_diagonalizeCImatrix(this) end do end if do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies + auxEnergy=0 + do sysI=1, this%numberOfDisplacedSystems + auxEnergy= auxEnergy+ & + this%configurationCoefficients%values(sysI,state)**2*& + this%configurationExchangeMatrix(speciesID)%values(sysI,sysI) + do sysII=sysI+1, this%numberOfDisplacedSystems + auxEnergy= auxEnergy+ & + 2.0_8*this%configurationCoefficients%values(sysI,state)*& + this%configurationCoefficients%values(sysII,state)*& + this%configurationExchangeMatrix(speciesID)%values(sysI,sysII) + end do + end do + write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(speciesID)%name ) // & + " Exchange energy = ", auxEnergy + do otherSpeciesID=speciesID, molecularSystem_instance%numberOfQuantumSpecies auxEnergy=0 do sysI=1, this%numberOfDisplacedSystems auxEnergy= auxEnergy+ & this%configurationCoefficients%values(sysI,state)**2*& - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(sysI,sysI) + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,sysI) do sysII=sysI+1, this%numberOfDisplacedSystems auxEnergy= auxEnergy+ & 2.0_8*this%configurationCoefficients%values(sysI,state)*& this%configurationCoefficients%values(sysII,state)*& - this%configurationTwoParticlesMatrix(speciesID,otherSpeciesID)%values(sysI,sysII) + this%configurationHartreeMatrix(speciesID,otherSpeciesID)%values(sysI,sysII) end do end do write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(speciesID)%name ) // & "/"//trim( MolecularSystem_instance%species(otherSpeciesID)%name ) // & - " Coulomb energy = ", auxEnergy + " Hartree energy = ", auxEnergy end do end do end do @@ -2121,86 +2186,76 @@ subroutine NonOrthogonalCI_diagonalizeCImatrix(this) end subroutine NonOrthogonalCI_diagonalizeCImatrix - subroutine NonOrthogonalCI_generateDensities(this) + !> + !! @brief Generates one molecular system combining all the displaced geometries and coefficients + !! + !! @param this + !< + subroutine NonOrthogonalCI_generateSuperposedSystem(this) implicit none type(NonOrthogonalCI) :: this type(MolecularSystem) :: auxMolecularSystem - type(Matrix), allocatable :: auxCoefficients(:),mergedCoefficients(:) - type(Matrix), allocatable :: overlapMatrix(:), kineticMatrix(:), attractionMatrix(:), externalPotMatrix(:) - type(IVector), allocatable :: sysBasisList(:,:),auxBasisList(:) - type(Matrix), allocatable :: mergedDensityMatrix(:,:) - - type(Matrix) :: auxMatrix, densityEigenVectors, auxdensityEigenVectors - type(Vector) :: densityEigenValues, auxdensityEigenValues + type(Matrix), allocatable :: auxCoefficients(:) + type(IVector), allocatable :: auxBasisList(:) - type(Matrix) :: molecularOverlapMatrix - type(Matrix), allocatable :: inverseOverlapMatrix(:) - integer :: state - integer :: i,ii,j,jj,k,mu,nu, sysI, sysII, speciesID - integer :: particlesPerOrbital + integer :: i, sysI, speciesID + integer :: numberOfSpecies - integer :: densUnit - character(100) :: densFile - character(50) :: arguments(2), auxString real(8) :: timeA - - !$ timeA = omp_get_wtime() + !$ timeA = omp_get_wtime() + if(CONTROL_instance%CI_STATES_TO_PRINT .eq. 0) return - allocate(mergedCoefficients(molecularSystem_instance%numberOfQuantumSpecies),& - auxCoefficients(molecularSystem_instance%numberOfQuantumSpecies),& - sysBasisList(this%numberOfDisplacedSystems,molecularSystem_instance%numberOfQuantumSpecies),& - auxBasisList(molecularSystem_instance%numberOfQuantumSpecies),& - overlapMatrix(molecularSystem_instance%numberOfQuantumSpecies),& - kineticMatrix(molecularSystem_instance%numberOfQuantumSpecies),& - attractionMatrix(molecularSystem_instance%numberOfQuantumSpecies),& - externalPotMatrix(molecularSystem_instance%numberOfQuantumSpecies),& - mergedDensityMatrix(CONTROL_instance%CI_STATES_TO_PRINT,molecularSystem_instance%numberOfQuantumSpecies)) + numberOfSpecies=molecularSystem_instance%numberOfQuantumSpecies + + allocate(this%sysBasisList(this%numberOfDisplacedSystems,numberOfSpecies),& + auxCoefficients(numberOfSpecies),& + auxBasisList(numberOfSpecies)) !Create a super molecular system !!!Merge coefficients from system 1 and system 2 - call MolecularSystem_mergeTwoSystems(molecularSystem_instance, this%MolecularSystems(1), this%MolecularSystems(2), & - sysBasisList(1,:),sysBasisList(2,:)) - + call MolecularSystem_mergeTwoSystems(this%mergedMolecularSystem, this%molecularSystems(1), this%molecularSystems(2), & + this%sysBasisList(1,:),this%sysBasisList(2,:)) + call NonOrthogonalCI_mergeCoefficients(this%HFCoefficients(1,:),this%HFCoefficients(2,:),& - this%MolecularSystems(1),this%MolecularSystems(2),molecularSystem_instance,& - sysBasisList(1,:),sysBasisList(2,:),mergedCoefficients) + this%molecularSystems(1),this%molecularSystems(2),this%mergedMolecularSystem,& + this%sysBasisList(1,:),this%sysBasisList(2,:),this%mergedCoefficients(:)) - ! do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - ! print *, "2", speciesID, "ocupationNumber", MolecularSystem_getOcupationNumber(speciesID) - ! print *, "2", speciesID, "mergedCoefficients" - ! call Matrix_show(mergedCoefficients(speciesID)) + ! do speciesID=1, numberOfSpecies + ! print *, "2", speciesID, "ocupationNumber", MolecularSystem_getOcupationNumber(speciesID,this%mergedMolecularSystem) + ! print *, "2", speciesID, "mergedCoefficients" + ! call Matrix_show(this%mergedCoefficients(speciesID)) ! end do ! !! Loop other systems expanding the merged coefficients matrix do sysI=3, this%numberOfDisplacedSystems - call MolecularSystem_copyConstructor(auxMolecularSystem,molecularSystem_instance) - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - call Matrix_copyConstructor(auxCoefficients(speciesID), mergedCoefficients(speciesID)) + call MolecularSystem_copyConstructor(auxMolecularSystem,this%mergedMolecularSystem) + do speciesID=1, numberOfSpecies + call Matrix_copyConstructor(auxCoefficients(speciesID), this%mergedCoefficients(speciesID)) end do - call MolecularSystem_mergeTwoSystems(molecularSystem_instance, auxMolecularSystem, this%MolecularSystems(sysI), & - auxBasisList,sysBasisList(sysI,:),reorder=.false.) + call MolecularSystem_mergeTwoSystems(this%mergedMolecularSystem, auxMolecularSystem, this%molecularSystems(sysI), & + auxBasisList,this%sysBasisList(sysI,:),reorder=.false.) call NonOrthogonalCI_mergeCoefficients(auxCoefficients,this%HFCoefficients(sysI,:),& - auxMolecularSystem,this%MolecularSystems(sysI),molecularSystem_instance,& - auxBasisList,sysBasisList(sysI,:),mergedCoefficients) - ! do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - ! print *, sysI, speciesID, "ocupationNumber", MolecularSystem_getOcupationNumber(speciesID) - ! print *, sysI, speciesID, "mergedCoefficients" - ! call Matrix_show(mergedCoefficients(speciesID)) + auxMolecularSystem,this%molecularSystems(sysI),this%mergedMolecularSystem,& + auxBasisList,this%sysBasisList(sysI,:),this%mergedCoefficients(:)) + ! do speciesID=1, numberOfSpecies + ! print *, sysI, speciesID, "ocupationNumber", MolecularSystem_getOcupationNumber(speciesID,this%mergedMolecularSystem) + ! print *, sysI, speciesID, "mergedCoefficients" + ! call Matrix_show(this%mergedCoefficients(speciesID)) ! end do end do !!!Fix basis list size do sysI=1, this%numberOfDisplacedSystems - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - call Vector_copyConstructorInteger(auxBasisList(speciesID),sysBasisList(sysI,speciesID)) - call Vector_constructorInteger(sysBasisList(sysI,speciesID), MolecularSystem_getTotalNumberOfContractions(speciesID), 0) + do speciesID=1, numberOfSpecies + call Vector_copyConstructorInteger(auxBasisList(speciesID),this%sysBasisList(sysI,speciesID)) + call Vector_constructorInteger(this%sysBasisList(sysI,speciesID), MolecularSystem_getTotalNumberOfContractions(speciesID,this%mergedMolecularSystem), 0) do i=1, size(auxBasisList(speciesID)%values) - sysBasisList(sysI,speciesID)%values(i)=auxBasisList(speciesID)%values(i) + this%sysBasisList(sysI,speciesID)%values(i)=auxBasisList(speciesID)%values(i) end do ! print *, "sysI", sysI, "speciesID", speciesID, "after list" - ! call Vector_showInteger(sysBasisList(sysI,speciesID)) + ! call Vector_showInteger(this%sysBasisList(sysI,speciesID)) end do end do @@ -2209,51 +2264,76 @@ subroutine NonOrthogonalCI_generateDensities(this) write(*,*) "---------------------------------- " ! call MolecularSystem_showInformation() ! call MolecularSystem_showParticlesInformation() + call MolecularSystem_copyConstructor(molecularSystem_instance,this%mergedMolecularSystem) call MolecularSystem_showCartesianMatrix(molecularSystem_instance) + particleManager_instance => molecularSystem_instance%allParticles + call ParticleManager_setOwner() call MolecularSystem_saveToFile() - - ! do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - ! write(*,*) "" - ! write(*,*) " Merged Occupied Eigenvectors for: ", trim( MolecularSystem_instance%species(speciesID)%name ) - ! write(*,*) "---------------------------------- " - ! write(*,*) "" - ! print *, "contractions", speciesID, int(MolecularSystem_getTotalNumberOfContractions(speciesID),8) - ! print *, "ocupation", speciesID, int(MolecularSystem_getOcupationNumber(speciesID),8) - ! call Matrix_constructor(auxCoefficients(speciesID),int(MolecularSystem_getTotalNumberOfContractions(speciesID),8),& - ! int(MolecularSystem_getOcupationNumber(speciesID),8),0.0_8) - ! do i=1, MolecularSystem_getTotalNumberOfContractions(speciesID) - ! do j=1, MolecularSystem_getOcupationNumber(speciesID) - ! auxCoefficients(speciesID)%values(i,j)=mergedCoefficients(speciesID)%values(i,j) - ! end do - ! end do - ! call Matrix_show(auxCoefficients(speciesID)) + + ! do speciesID=1, numberOfSpecies + ! write(*,*) "" + ! write(*,*) " Merged Occupied Eigenvectors for: ", trim( MolecularSystem_instance%species(speciesID)%name ) + ! write(*,*) "---------------------------------- " + ! write(*,*) "" + ! print *, "contractions", speciesID, int(MolecularSystem_getTotalNumberOfContractions(speciesID),8) + ! print *, "ocupation", speciesID, int(MolecularSystem_getOcupationNumber(speciesID),8) + ! call Matrix_constructor(auxCoefficients(speciesID),int(MolecularSystem_getTotalNumberOfContractions(speciesID),8),& + ! int(MolecularSystem_getOcupationNumber(speciesID),8),0.0_8) + ! do i=1, MolecularSystem_getTotalNumberOfContractions(speciesID) + ! do j=1, MolecularSystem_getOcupationNumber(speciesID) + ! auxCoefficients(speciesID)%values(i,j)=mergedCoefficients(speciesID)%values(i,j) + ! end do + ! end do + ! call Matrix_show(auxCoefficients(speciesID)) ! end do !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time creating supermolecular system : ", omp_get_wtime() - timeA ," (s)" !$ timeA = omp_get_wtime() + + deallocate(auxCoefficients,& + auxBasisList) + + return - print *, "" - print *, "Computing overlap integrals in for the superposed systems..." - print *, "" - !!Compute overlap integrals - do speciesID = 1, MolecularSystem_instance%numberOfQuantumSpecies - call DirectIntegralManager_getOverlapIntegrals(molecularSystem_instance,speciesID,overlapMatrix(speciesID)) - ! call DirectIntegralManager_getKineticIntegrals(molecularSystem_instance,speciesID,kineticMatrix(speciesID)) - ! if ( CONTROL_instance%REMOVE_TRANSLATIONAL_CONTAMINATION ) then - ! kineticMatrix(speciesID)%values = & - ! kineticMatrix(speciesID)%values * & - ! ( 1.0_8/MolecularSystem_getMass( speciesID ) -1.0_8 / MolecularSystem_getTotalMass() ) - ! else - ! kineticMatrix(speciesID)%values = & - ! kineticMatrix(speciesID)%values / & - ! MolecularSystem_getMass( speciesID ) - ! end if + end subroutine NonOrthogonalCI_generateSuperposedSystem - ! call DirectIntegralManager_getAttractionIntegrals(molecularSystem_instance,speciesID,attractionMatrix(speciesID)) - ! attractionMatrix(speciesID)%values=attractionMatrix(speciesID)%values*(-MolecularSystem_getCharge(speciesID)) + !> + !! @brief Generates the NOCI density matrix in the superposed molecular system + !! + !! @param this + !< + subroutine NonOrthogonalCI_buildDensityMatrix(this) + implicit none + type(NonOrthogonalCI) :: this + + type(Matrix) :: molecularOverlapMatrix + type(Matrix), allocatable :: inverseOverlapMatrix(:) !,kineticMatrix(:), attractionMatrix(:), externalPotMatrix(:) + integer :: state + integer :: i,ii,j,jj,mu,nu, sysI, sysII, speciesID, otherSpeciesID + integer :: particlesPerOrbital + integer :: numberOfSpecies + + integer :: densUnit + character(100) :: densFile + character(50) :: arguments(2), auxString + type(Matrix), allocatable :: exchangeCorrelationMatrices(:) + type(Matrix) :: dftEnergyMatrix + real(8), allocatable :: particlesInGrid(:) + real(8) :: timeA - ! if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & - ! call DirectIntegralManager_getExternalPotentialIntegrals(molecularSystem_instance,speciesID,externalPotMatrix(speciesID)) + !$ timeA = omp_get_wtime() + + if(CONTROL_instance%CI_STATES_TO_PRINT .eq. 0) return + + numberOfSpecies=molecularSystem_instance%numberOfQuantumSpecies + + allocate(InverseOverlapMatrix(numberOfSpecies)) + + print *, "" + print *, "Computing overlap integrals for the superposed systems..." + print *, "" + do speciesID = 1, numberOfSpecies + call DirectIntegralManager_getOverlapIntegrals(molecularSystem_instance,speciesID,this%mergedOverlapMatrix(speciesID)) end do !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for supermolecular 1-body integrals : ", omp_get_wtime() - timeA ," (s)" !$ timeA = omp_get_wtime() @@ -2261,30 +2341,29 @@ subroutine NonOrthogonalCI_generateDensities(this) print *, "" print *, "Building merged density matrices for the superposed systems..." print *, "" - allocate(InverseOverlapMatrix(molecularSystem_instance%numberOfQuantumSpecies)) !!Build the merged density matrix do state=1, CONTROL_instance%CI_STATES_TO_PRINT - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - call Matrix_constructor(mergedDensityMatrix(state,speciesID), int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & + do speciesID=1, numberOfSpecies + call Matrix_constructor(this%mergedDensityMatrix(state,speciesID), int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), 0.0_8) end do end do !!Fill the merged density matrix ! "Diagonal" terms - same system do sysI=1, this%numberOfDisplacedSystems - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - particlesPerOrbital=MolecularSystem_getEta(speciesID,this%MolecularSystems(sysI)) + do speciesID=1, numberOfSpecies + particlesPerOrbital=MolecularSystem_getEta(speciesID,this%molecularSystems(sysI)) do mu = 1 , MolecularSystem_getTotalNumberOfContractions(speciesID) - if(sysBasisList(sysI,speciesID)%values(mu) .eq. 0) cycle + if(this%sysBasisList(sysI,speciesID)%values(mu) .eq. 0) cycle do nu = 1 , MolecularSystem_getTotalNumberOfContractions(speciesID) - if(sysBasisList(sysI,speciesID)%values(nu) .eq. 0) cycle - do i = 1 , MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) - ii=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI))*(sysI-1)+i + if(this%sysBasisList(sysI,speciesID)%values(nu) .eq. 0) cycle + do i = 1 , MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) + ii=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI))*(sysI-1)+i do state=1, CONTROL_instance%CI_STATES_TO_PRINT - mergedDensityMatrix(state,speciesID)%values(mu,nu) = mergedDensityMatrix(state,speciesID)%values(mu,nu) + & + this%mergedDensityMatrix(state,speciesID)%values(mu,nu) = this%mergedDensityMatrix(state,speciesID)%values(mu,nu) + & this%configurationCoefficients%values(sysI,state)**2*& - mergedCoefficients(speciesID)%values(mu,ii)*& - mergedCoefficients(speciesID)%values(nu,ii)*& + this%mergedCoefficients(speciesID)%values(mu,ii)*& + this%mergedCoefficients(speciesID)%values(nu,ii)*& particlesPerOrbital end do end do @@ -2297,64 +2376,64 @@ subroutine NonOrthogonalCI_generateDensities(this) do sysII=sysI+1, this%numberOfDisplacedSystems if( abs(this%configurationOverlapMatrix%values(sysI,sysII)) .lt. CONTROL_instance%CONFIGURATION_OVERLAP_THRESHOLD ) cycle !!Compute molecular overlap matrix and its inverse - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies + do speciesID=1, numberOfSpecies call Matrix_constructor(molecularOverlapMatrix, & - int(MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)),8), & - int(MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysII)),8), 0.0_8 ) + int(MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)),8), & + int(MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysII)),8), 0.0_8 ) call Matrix_constructor(inverseOverlapMatrix(speciesID), & - int(MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)),8), & - int(MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysII)),8), 0.0_8 ) + int(MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)),8), & + int(MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysII)),8), 0.0_8 ) do mu=1, MolecularSystem_getTotalNumberOfContractions(speciesID) !sysI - if(sysBasisList(sysI,speciesID)%values(mu) .eq. 0) cycle + if(this%sysBasisList(sysI,speciesID)%values(mu) .eq. 0) cycle do nu=1, MolecularSystem_getTotalNumberOfContractions(speciesID) !sysII - if(sysBasisList(sysII,speciesID)%values(nu) .eq. 0) cycle - do i = 1 , MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) - ii=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI))*(sysI-1)+i - do j = 1 , MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysII)) - jj=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysII))*(sysII-1)+j + if(this%sysBasisList(sysII,speciesID)%values(nu) .eq. 0) cycle + do i = 1 , MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) + ii=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI))*(sysI-1)+i + do j = 1 , MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysII)) + jj=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysII))*(sysII-1)+j ! print *, "i, j, mu, nu, coefI, coefII", i,j,mu,nu,mergedCoefficients(speciesID)%values(mu,ii),mergedCoefficients(speciesID)%values(nu,jj) molecularOverlapMatrix%values(i,j)=molecularOverlapMatrix%values(i,j)+& - mergedCoefficients(speciesID)%values(mu,ii)*& - mergedCoefficients(speciesID)%values(nu,jj)*& - overlapMatrix(speciesID)%values(mu,nu) + this%mergedCoefficients(speciesID)%values(mu,ii)*& + this%mergedCoefficients(speciesID)%values(nu,jj)*& + this%mergedOverlapMatrix(speciesID)%values(mu,nu) end do end do end do end do ! print *, "molecularOverlapMatrix sysI, sysII, speciesID", sysI, sysII, speciesID ! call Matrix_show(molecularOverlapMatrix) - if(MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) .ne. 0) & + if(MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) .ne. 0) & inverseOverlapMatrix(speciesID)=Matrix_inverse(molecularOverlapMatrix) end do ! Compute density contributions - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - particlesPerOrbital=MolecularSystem_getEta(speciesID,this%MolecularSystems(sysI)) + do speciesID=1, numberOfSpecies + particlesPerOrbital=MolecularSystem_getEta(speciesID,this%molecularSystems(sysI)) do mu = 1 , MolecularSystem_getTotalNumberOfContractions(speciesID) - if(sysBasisList(sysI,speciesID)%values(mu) .eq. 0) cycle + if(this%sysBasisList(sysI,speciesID)%values(mu) .eq. 0) cycle do nu = 1 , MolecularSystem_getTotalNumberOfContractions(speciesID) - if(sysBasisList(sysII,speciesID)%values(nu) .eq. 0) cycle + if(this%sysBasisList(sysII,speciesID)%values(nu) .eq. 0) cycle do state=1, CONTROL_instance%CI_STATES_TO_PRINT - do i = 1 , MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI)) - ii=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysI))*(sysI-1)+i - do j = 1 , MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysII)) - jj=MolecularSystem_getOcupationNumber(speciesID,this%MolecularSystems(sysII))*(sysII-1)+j - mergedDensityMatrix(state,speciesID)%values(mu,nu) = mergedDensityMatrix(state,speciesID)%values(mu,nu) + & + do i = 1 , MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) + ii=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI))*(sysI-1)+i + do j = 1 , MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysII)) + jj=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysII))*(sysII-1)+j + this%mergedDensityMatrix(state,speciesID)%values(mu,nu) = this%mergedDensityMatrix(state,speciesID)%values(mu,nu) + & this%configurationCoefficients%values(sysI,state)*& this%configurationCoefficients%values(sysII,state)*& this%configurationOverlapMatrix%values(sysI,sysII)*& inverseOverlapMatrix(speciesID)%values(j,i)*& - mergedCoefficients(speciesID)%values(mu,ii)*& - mergedCoefficients(speciesID)%values(nu,jj)*& + this%mergedCoefficients(speciesID)%values(mu,ii)*& + this%mergedCoefficients(speciesID)%values(nu,jj)*& particlesPerOrbital - mergedDensityMatrix(state,speciesID)%values(nu,mu) = mergedDensityMatrix(state,speciesID)%values(nu,mu) + & + this%mergedDensityMatrix(state,speciesID)%values(nu,mu) = this%mergedDensityMatrix(state,speciesID)%values(nu,mu) + & this%configurationCoefficients%values(sysI,state)*& this%configurationCoefficients%values(sysII,state)*& this%configurationOverlapMatrix%values(sysI,sysII)*& inverseOverlapMatrix(speciesID)%values(j,i)*& - mergedCoefficients(speciesID)%values(mu,ii)*& - mergedCoefficients(speciesID)%values(nu,jj)*& + this%mergedCoefficients(speciesID)%values(mu,ii)*& + this%mergedCoefficients(speciesID)%values(nu,jj)*& particlesPerOrbital end do end do @@ -2365,7 +2444,7 @@ subroutine NonOrthogonalCI_generateDensities(this) !!symmetrize ! do mu = 1 , MolecularSystem_getTotalNumberOfContractions(speciesID) ! do nu = mu+1 , MolecularSystem_getTotalNumberOfContractions(speciesID) - ! mergedDensityMatrix(state,speciesID)%values(nu,mu) = mergedDensityMatrix(state,speciesID)%values(mu,nu) + ! this%mergedDensityMatrix(state,speciesID)%values(nu,mu) = this%mergedDensityMatrix(state,speciesID)%values(mu,nu) ! end do ! end do end do @@ -2377,166 +2456,245 @@ subroutine NonOrthogonalCI_generateDensities(this) densFile = trim(CONTROL_instance%INPUT_FILE)//"Matrices.ci" open(unit = densUnit, file=trim(densFile), status="replace", form="formatted") do state=1, CONTROL_instance%CI_STATES_TO_PRINT - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies - ! print *, "mergedDensityMatrix", state, trim( MolecularSystem_instance%species(speciesID)%name ) - ! call Matrix_show(mergedDensityMatrix(state,speciesID)) + do speciesID=1, numberOfSpecies + ! print *, "this%mergedDensityMatrix", state, trim( MolecularSystem_instance%species(speciesID)%name ) + ! call Matrix_show(this%mergedDensityMatrix(state,speciesID)) write(auxString,*) state arguments(2) = trim(MolecularSystem_instance%species(speciesID)%name) arguments(1) = "DENSITYMATRIX"//trim(adjustl(auxString)) - call Matrix_writeToFile ( mergedDensityMatrix(state,speciesID), densUnit , arguments=arguments(1:2) ) + call Matrix_writeToFile ( this%mergedDensityMatrix(state,speciesID), densUnit , arguments=arguments(1:2) ) end do end do + if(CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL.ne."NONE" .or. & + CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL.ne."NONE") then + print *, "Superposed DFT energies:" + + allocate(exchangeCorrelationMatrices(numberOfSpecies), & + particlesInGrid(numberOfSpecies)) + call DensityFunctionalTheory_buildFinalGrid() + do state=1, CONTROL_instance%CI_STATES_TO_PRINT + call Matrix_constructor(dftEnergyMatrix, int(numberOfSpecies,8), & + int(numberOfSpecies,8), 0.0_8 ) + do speciesID=1, numberOfSpecies + call Matrix_constructor(exchangeCorrelationMatrices(speciesID), int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), 0.0_8) + end do + call DensityFunctionalTheory_finalDFT(this%mergedDensityMatrix(state,1:numberOfSpecies), & + exchangeCorrelationMatrices, & + dftEnergyMatrix, & + particlesInGrid) + + do speciesID = 1, MolecularSystem_instance%numberOfQuantumSpecies + do otherSpeciesID = speciesID, MolecularSystem_instance%numberOfQuantumSpecies + write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(speciesID)%name ) // & + "/"//trim( MolecularSystem_instance%species(otherSpeciesID)%name ) // & + " DFT Corr. energy = ", dftEnergyMatrix%values(speciesID,otherSpeciesID) + end do + end do + end do + end if + + close(densUnit) + + deallocate(inverseOverlapMatrix) + !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for merging density matrices : ", omp_get_wtime() - timeA ," (s)" - !$ timeA = omp_get_wtime() + + return + ! allocate(kineticMatrix(numberOfSpecies),& + ! attractionMatrix(numberOfSpecies),& + ! externalPotMatrix(numberOfSpecies)) + ! do speciesID = 1, numberOfSpecies + ! call DirectIntegralManager_getKineticIntegrals(molecularSystem_instance,speciesID,kineticMatrix(speciesID)) + ! if ( CONTROL_instance%REMOVE_TRANSLATIONAL_CONTAMINATION ) then + ! kineticMatrix(speciesID)%values = & + ! kineticMatrix(speciesID)%values * & + ! ( 1.0_8/MolecularSystem_getMass( speciesID ) -1.0_8 / MolecularSystem_getTotalMass() ) + ! else + ! kineticMatrix(speciesID)%values = & + ! kineticMatrix(speciesID)%values / & + ! MolecularSystem_getMass( speciesID ) + ! end if + + ! call DirectIntegralManager_getAttractionIntegrals(molecularSystem_instance,speciesID,attractionMatrix(speciesID)) + ! attractionMatrix(speciesID)%values=attractionMatrix(speciesID)%values*(-MolecularSystem_getCharge(speciesID)) + + ! if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & + ! call DirectIntegralManager_getExternalPotentialIntegrals(molecularSystem_instance,speciesID,externalPotMatrix(speciesID)) + ! end do ! write(*,*) "" ! write(*,*) "==========================================================" ! write(*,*) " ONE BODY ENERGY CONTRIBUTIONS OF THE SUPERPOSED SYSTEMS: " ! write(*,*) "" ! do state=1, CONTROL_instance%CI_STATES_TO_PRINT ! write(*,*) " STATE: ", state - ! do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies + ! do speciesID=1, numberOfSpecies ! write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(speciesID)%name ) // & - ! " Kinetic energy = ", sum(transpose(mergedDensityMatrix(state,speciesID)%values)*kineticMatrix(speciesID)%values) + ! " Kinetic energy = ", sum(transpose(this%mergedDensityMatrix(state,speciesID)%values)*kineticMatrix(speciesID)%values) ! write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(speciesID)%name ) // & - ! "/Fixed interact. energy = ", sum(transpose(mergedDensityMatrix(state,speciesID)%values)*attractionMatrix(speciesID)%values) + ! "/Fixed interact. energy = ", sum(transpose(this%mergedDensityMatrix(state,speciesID)%values)*attractionMatrix(speciesID)%values) ! if( CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & ! write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(speciesID)%name) // & - ! " Ext Pot energy = ", sum(transpose(mergedDensityMatrix(state,speciesID)%values)*externalPotMatrix(speciesID)%values) + ! " Ext Pot energy = ", sum(transpose(this%mergedDensityMatrix(state,speciesID)%values)*externalPotMatrix(speciesID)%values) ! print *, "" ! end do ! print *, "" ! end do + ! deallocate(kineticMatrix,& + ! attractionMatrix,& + ! externalPotMatrix) - !! Natural orbitals + end subroutine NonOrthogonalCI_buildDensityMatrix - if (CONTROL_instance%CI_NATURAL_ORBITALS) then - write(*,*) "" - write(*,*) "=============================================" - write(*,*) " NATURAL ORBITALS OF THE SUPERPOSED SYSTEMS: " - write(*,*) "" + !> + !! @brief Generates the NOCI natural orbitals from the NOCI density matrix in the superposed molecular system + !! + !! @param this + !< + subroutine NonOrthogonalCI_getNaturalOrbitals(this) + implicit none + type(NonOrthogonalCI) :: this - do state=1, CONTROL_instance%CI_STATES_TO_PRINT + type(Matrix) :: auxMatrix, densityEigenVectors, auxdensityEigenVectors + type(Vector) :: auxVector, densityEigenValues, auxdensityEigenValues - write(*,*) " STATE: ", state + integer :: state + integer :: i,j,k,speciesID + integer :: numberOfSpecies - do speciesID=1, molecularSystem_instance%numberOfQuantumSpecies + integer :: densUnit + character(100) :: densFile + character(50) :: arguments(2), auxString + real(8) :: timeA - write(*,*) "" - write(*,*) " Natural Orbitals in state: ", state, " for: ", trim( MolecularSystem_instance%species(speciesID)%name ) - write(*,*) "-----------------" + !$ timeA = omp_get_wtime() + if(.not. CONTROL_instance%CI_NATURAL_ORBITALS) return + if(CONTROL_instance%CI_STATES_TO_PRINT .eq. 0) return - call Vector_constructor ( densityEigenValues, & - int(MolecularSystem_getTotalNumberOfContractions(speciesID),4), 0.0_8 ) - call Matrix_constructor ( densityEigenVectors, & - int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & - int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), 0.0_8 ) + numberOfSpecies=molecularSystem_instance%numberOfQuantumSpecies - call Vector_constructor ( auxdensityEigenValues, & - int(MolecularSystem_getTotalNumberOfContractions(speciesID),4), 0.0_8 ) - call Matrix_constructor ( auxdensityEigenVectors, & - int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & - int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), 0.0_8 ) + write(*,*) "" + write(*,*) "=============================================" + write(*,*) " NATURAL ORBITALS OF THE SUPERPOSED SYSTEMS: " + write(*,*) "" + !! Open file - to write density matrices + densUnit = 29 - ! print *,"Matriz de overlap " - ! call Matrix_show( overlapMatrix(speciesID) ) + densFile = trim(CONTROL_instance%INPUT_FILE)//"Matrices.ci" + open(unit = densUnit, file=trim(densFile), status="old", form="formatted", position="append") - !! Lowdin orthogonalization of the density matrix - auxMatrix = Matrix_pow( overlapMatrix(speciesID), 0.5_8 ) + do state=1, CONTROL_instance%CI_STATES_TO_PRINT - auxMatrix%values=matmul(matmul(auxMatrix%values,mergedDensityMatrix(state,speciesID)%values),auxMatrix%values) + write(*,*) " STATE: ", state - ! print *, "Diagonalizing non orthogonal CI density Matrix..." + do speciesID=1, numberOfSpecies - !! Calcula valores y vectores propios de matriz de densidad CI ortogonal. - call Matrix_eigen(auxMatrix , auxdensityEigenValues, auxdensityEigenVectors, SYMMETRIC ) + write(*,*) "" + write(*,*) " Natural Orbitals in state: ", state, " for: ", trim( MolecularSystem_instance%species(speciesID)%name ) + write(*,*) "--------------------------------------------------------------" - !! Transform back to the atomic basis - auxMatrix = Matrix_pow( overlapMatrix(speciesID), -0.5_8 ) + call Vector_constructor ( densityEigenValues, & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),4), 0.0_8 ) + call Matrix_constructor ( densityEigenVectors, & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), 0.0_8 ) - auxdensityEigenVectors%values=matmul(auxMatrix%values,auxdensityEigenVectors%values) + call Vector_constructor ( auxdensityEigenValues, & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),4), 0.0_8 ) + call Matrix_constructor ( auxdensityEigenVectors, & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), & + int(MolecularSystem_getTotalNumberOfContractions(speciesID),8), 0.0_8 ) - ! reorder and count significant occupations - k=0 - do i = 1, MolecularSystem_getTotalNumberOfContractions(speciesID) - densityEigenValues%values(i) = auxdensityEigenValues%values(MolecularSystem_getTotalNumberOfContractions(speciesID) - i + 1) - densityEigenVectors%values(:,i) = auxdensityEigenVectors%values(:,MolecularSystem_getTotalNumberOfContractions(speciesID) - i + 1) - if(densityEigenValues%values(i) .ge. 0.0000000000001 ) k=k+1 - end do - if(k .eq. 0) k=1 - ! Print eigenvectors with occupation larger than 0.01 - call Matrix_constructor(auxMatrix,int(MolecularSystem_getTotalNumberOfContractions(speciesID),8),int(k,8),0.0_8) - do i=1, MolecularSystem_getTotalNumberOfContractions(speciesID) - do j=1, k - auxMatrix%values(i,j)=densityEigenVectors%values(i,j) + !! Lowdin orthogonalization of the density matrix + auxMatrix = Matrix_pow(this%mergedOverlapMatrix(speciesID), 0.5_8, method="SVD" ) + + auxMatrix%values=matmul(matmul(auxMatrix%values,this%mergedDensityMatrix(state,speciesID)%values),auxMatrix%values) + + ! print *, "Diagonalizing non orthogonal CI density Matrix..." + + !! Calcula valores y vectores propios de matriz de densidad CI ortogonal. + call Matrix_eigen(auxMatrix , auxdensityEigenValues, auxdensityEigenVectors, SYMMETRIC ) + + !! Transform back to the atomic basis + auxMatrix = Matrix_pow(this%mergedOverlapMatrix(speciesID), -0.5_8, method="SVD" ) + + auxdensityEigenVectors%values=matmul(auxMatrix%values,auxdensityEigenVectors%values) + + ! reorder and count significant occupations + k=0 + do i = 1, MolecularSystem_getTotalNumberOfContractions(speciesID) + densityEigenValues%values(i) = auxdensityEigenValues%values(MolecularSystem_getTotalNumberOfContractions(speciesID) - i + 1) + densityEigenVectors%values(:,i) = auxdensityEigenVectors%values(:,MolecularSystem_getTotalNumberOfContractions(speciesID) - i + 1) + if(abs(densityEigenValues%values(i)) .ge. 1.0E-4_8 ) k=k+1 + end do + if(k .eq. 0) k=1 + ! Print eigenvectors with occupation larger than 0.01 + call Vector_constructor(auxVector,k,0.0_8) + call Matrix_constructor(auxMatrix,int(MolecularSystem_getTotalNumberOfContractions(speciesID),8),int(k,8),0.0_8) + k=0 + do i=1, MolecularSystem_getTotalNumberOfContractions(speciesID) + if(abs(densityEigenValues%values(i)) .ge. 1.0E-4_8 ) then + k=k+1 + auxVector%values(k)=densityEigenValues%values(i) + do j=1, MolecularSystem_getTotalNumberOfContractions(speciesID) + auxMatrix%values(j,k)=densityEigenVectors%values(j,i) end do - end do - !densityEigenVectors - call Matrix_show( auxMatrix , & - rowkeys = MolecularSystem_getlabelsofcontractions( speciesID ), & - columnkeys = string_convertvectorofrealstostring( densityEigenValues ),& - flags=WITH_BOTH_KEYS) - - write(*,"(A10,A10,A20,I5,A15,F17.12)") "number of ", trim(MolecularSystem_getNameOfSpecie( speciesID )) ," particles in state", state , & - " density matrix: ", sum( transpose(mergedDensityMatrix(state,speciesID)%values)*overlapMatrix(speciesID)%values) - write(*,"(A10,A10,A40,F17.12)") "sum of ", trim(MolecularSystem_getNameOfSpecie( speciesID )) , "natural orbital occupations", sum(densityEigenValues%values) - - ! density matrix check - ! auxMatrix%values=0.0 - ! do mu=1, MolecularSystem_getTotalNumberOfContractions(speciesID) - ! do nu=1, MolecularSystem_getTotalNumberOfContractions(speciesID) - ! do k=1, MolecularSystem_getTotalNumberOfContractions(speciesID) - ! auxMatrix%values(mu,nu)=auxMatrix%values(mu,nu)+densityEigenValues%values(k)*& - ! densityEigenVectors%values(mu,k)*densityEigenVectors%values(nu,k) - ! end do - ! end do - ! end do - ! print *, "atomicDensityMatrix again" - ! call Matrix_show(auxMatrix) - - write(auxString,*) state - - arguments(2) = trim( MolecularSystem_instance%species(speciesID)%name ) - arguments(1) = "NATURALORBITALS"//trim(adjustl(auxstring)) - - call Matrix_writeToFile ( densityEigenVectors, densUnit , arguments=arguments(1:2) ) - - arguments(2) = trim( MolecularSystem_instance%species(speciesID)%name ) - arguments(1) = "OCCUPATIONS"//trim(adjustl(auxstring)) - - call Vector_writeToFile( densityEigenValues, densUnit, arguments=arguments(1:2) ) - - write(*,*) " End of natural orbitals in state: ", state, " for: ", trim( MolecularSystem_instance%species(speciesID)%name ) + end if end do + !densityEigenVectors + call Matrix_show( auxMatrix , & + rowkeys = MolecularSystem_getlabelsofcontractions( speciesID ), & + columnkeys = string_convertvectorofrealstostring( auxVector ),& + flags=WITH_BOTH_KEYS) + + write(*,"(A10,A10,A20,I5,A15,F17.12)") "number of ", trim(MolecularSystem_getNameOfSpecie( speciesID )) ," particles in state", state , & + " density matrix: ", sum( transpose(this%mergedDensityMatrix(state,speciesID)%values)*this%mergedOverlapMatrix(speciesID)%values) + write(*,"(A10,A10,A40,F17.12)") "sum of ", trim(MolecularSystem_getNameOfSpecie( speciesID )) , "natural orbital occupations", sum(densityEigenValues%values) + + ! density matrix check + ! auxMatrix%values=0.0 + ! do mu=1, MolecularSystem_getTotalNumberOfContractions(speciesID) + ! do nu=1, MolecularSystem_getTotalNumberOfContractions(speciesID) + ! do k=1, MolecularSystem_getTotalNumberOfContractions(speciesID) + ! auxMatrix%values(mu,nu)=auxMatrix%values(mu,nu)+densityEigenValues%values(k)*& + ! densityEigenVectors%values(mu,k)*densityEigenVectors%values(nu,k) + ! end do + ! end do + ! end do + ! print *, "atomicDensityMatrix again" + ! call Matrix_show(auxMatrix) + + write(auxString,*) state + + arguments(2) = trim( MolecularSystem_instance%species(speciesID)%name ) + arguments(1) = "NATURALORBITALS"//trim(adjustl(auxstring)) + + call Matrix_writeToFile ( densityEigenVectors, densUnit , arguments=arguments(1:2) ) + + arguments(2) = trim( MolecularSystem_instance%species(speciesID)%name ) + arguments(1) = "OCCUPATIONS"//trim(adjustl(auxstring)) + + call Vector_writeToFile( densityEigenValues, densUnit, arguments=arguments(1:2) ) + + write(*,*) " End of natural orbitals in state: ", state, " for: ", trim( MolecularSystem_instance%species(speciesID)%name ) end do + end do + + write(*,*) "" + write(*,*) " END OF NATURAL ORBITALS" + write(*,*) "==============================" + write(*,*) "" - write(*,*) "" - write(*,*) " END OF NATURAL ORBITALS" - write(*,*) "==============================" - write(*,*) "" - end if - close(densUnit) - !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for NOCI natural orbitals : ", omp_get_wtime() - timeA ," (s)" - - deallocate(mergedCoefficients,& - auxCoefficients,& - sysBasisList,& - auxBasisList,& - overlapMatrix,& - kineticMatrix,& - attractionMatrix,& - externalPotMatrix,& - inverseOverlapMatrix,& - mergedDensityMatrix) - + return - - end subroutine NonOrthogonalCI_generateDensities + end subroutine NonOrthogonalCI_getNaturalOrbitals + !> !! @brief Calculate and Transform the four center integrals in one sweep without writing anything to disk !! @@ -2783,6 +2941,510 @@ subroutine NonOrthogonalCI_transformIntegralsMemory(mergedMolecularSystem, merge end subroutine NonOrthogonalCI_transformIntegralsMemory + + !> + !! @brief Save NOCI results to file + !! + !! @param + !< + subroutine NonOrthogonalCI_saveToFile(this) + type(NonOrthogonalCI) :: this + integer :: nociUnit, speciesID, numberOfSpecies, sysI + character(100) :: prefix, nociFile + character(50) :: arguments(2), auxString + + !Save merged molecular system + call MolecularSystem_copyConstructor(molecularSystem_instance,this%mergedMolecularSystem) + + prefix=trim(CONTROL_instance%INPUT_FILE)//"NOCI" + call MolecularSystem_saveToFile(prefix) + + numberOfSpecies=molecularSystem_instance%numberOfQuantumSpecies + + nociUnit=123 + nociFile=trim(prefix)//".states" + open(unit = nociUnit, file=trim(nociFile), status="replace", form="unformatted") + + arguments(1:1) = "NOCI-NUMBEROFDISPLACEDSYSTEMS" + call Vector_writeToFileInteger(unit=nociUnit, binary=.true., value=this%numberOfDisplacedSystems, arguments=arguments(1:1) ) + + arguments(1:1) = "NOCI-NUMBEROFSPECIES" + call Vector_writeToFileInteger(unit=nociUnit, binary=.true., value=numberOfSpecies, arguments=arguments(1:1) ) + + arguments(1:1) = "NOCI-CONFIGURATIONCOEFFICIENTS" + call Matrix_writeToFile ( this%configurationCoefficients, nociUnit , binary=.true., arguments=arguments(1:1) ) + + arguments(1:1) = "NOCI-CONFIGURATIONENERGIES" + call Vector_writeToFile ( this%statesEigenvalues, nociUnit , binary=.true., arguments=arguments(1:1) ) + + arguments(1) = "MERGEDCOEFFICIENTS" + do speciesID=1, numberOfSpecies + arguments(2) = trim(MolecularSystem_instance%species(speciesID)%name) + call Matrix_writeToFile ( this%mergedCoefficients(speciesID), nociUnit, binary=.true. , arguments=arguments(1:2) ) + end do + + do sysI=1, this%numberOfDisplacedSystems + do speciesID=1, numberOfSpecies + write(auxString,*) sysI + arguments(1) = "SYSBASISLIST"//trim(auxString) + arguments(2) = trim(MolecularSystem_instance%species(speciesID)%name) + call Vector_writeToFileInteger(this%sysBasisList(sysI,speciesID), nociUnit, binary=.true., arguments=arguments(1:2) ) + end do + end do + + ! do state=1, min(CONTROL_instance%NUMBER_OF_CI_STATES,this%numberOfDisplacedSystems) + ! end do + + ! do state=1, CONTROL_instance%CI_STATES_TO_PRINT + ! write(auxString,*) state + ! call Matrix_writeToFile ( this%mergedDensityMatrix(state,speciesID), densUnit , arguments=arguments(1:2) ) + ! end do + ! end do + + close(nociUnit) + + end subroutine NonOrthogonalCI_saveToFile + + !> + !! @brief Compute Franck-Condon factors from the current NOCI calculations and previous results read from file + !! + !! @param + !< + subroutine NonOrthogonalCI_computeFranckCondon(this) + type(NonOrthogonalCI) :: this + integer :: nociUnit, numberOfSpecies, occupationNumber,numberOfDisplacedSystems, numberOfContractions, dim2 + character(100) :: nociFile + type(Matrix) :: ciCoefficients + type(Vector) :: ciEnergies + type(Matrix), allocatable :: auxCoefficients(:), superMergedCoefficients(:) + type(IVector), allocatable :: sysListCur(:,:), sysListRef(:,:), orbListI(:), orbListII(:) + type(IVector) :: auxIVector + type(MolecularSystem) :: superMergedMolecularSystem + logical :: existFile + type(Matrix) :: molecularOverlapMatrix + type(Matrix), allocatable :: superOverlapMatrix(:), superMomentMatrix(:,:), inverseOverlapMatrix(:), molecularMomentMatrix(:,:) !,attractionMatrix(:), externalPotMatrix(:) + integer :: stateI, stateII + integer :: i,ii,j,jj,k,mu,nu,mumu,nunu,sysI, sysII, speciesID, otherSpeciesID + integer :: particlesPerOrbital + real(8) :: overlapDeterminant, trololo, trolololo(3), pointchargesdipole(3) + + integer :: densUnit + character(100) :: densFile + character(50) :: arguments(2), auxString + type(Matrix), allocatable :: franckCondonMatrix(:), transitionDipoleMatrix(:,:), refCurOverlapMatrix(:), refCurMomentMatrix(:,:) + type(Matrix) :: refCurTotalOverlap + real(8) :: timeA + + !$ timeA = omp_get_wtime() + + existFile=.false. + + nociFile = trim(CONTROL_instance%INPUT_FILE)//"refNOCI" + inquire( FILE = trim(nociFile)//".sys", EXIST = existFile ) + + if(.not. existFile) return + print *, "Found a reference molecular system for NOCI calculations ", trim(nociFile)//".sys" + + pointchargesdipole=0.0 + do i=1, size( MolecularSystem_instance%pointCharges ) + pointchargesdipole = pointchargesdipole + MolecularSystem_instance%pointCharges(i)%origin(:) * MolecularSystem_instance%pointCharges(i)%charge + end do + + + call MolecularSystem_loadFromFile("LOWDIN.SYS",nociFile) + call MolecularSystem_showInformation() + call MolecularSystem_showParticlesInformation() + call MolecularSystem_showCartesianMatrix() + + nociFile = trim(CONTROL_instance%INPUT_FILE)//"refNOCI.states" + inquire( FILE = trim(nociFile), EXIST = existFile ) + + if(.not. existFile) then + print *, "Did not find reference states for NOCI calculations ", nociFile + return + end if + print *, "Found reference states for NOCI calculations ", nociFile + print *, "Computing the Franck-Condon factors with respect to that system" + + nociUnit=123 + open(unit = nociUnit, file=trim(nociFile), status="old", form="unformatted") + + arguments(1) = "NOCI-NUMBEROFSPECIES" + call Vector_getFromFileInteger(1,unit=nociUnit, binary=.true., value=numberOfSpecies, arguments=arguments(1:1) ) + + arguments(1) = "NOCI-NUMBEROFDISPLACEDSYSTEMS" + call Vector_getFromFileInteger(1,unit=nociUnit, binary=.true., value=numberOfDisplacedSystems, arguments=arguments(1:1) ) + + allocate(auxCoefficients(numberOfSpecies)) + allocate(sysListCur(numberOfDisplacedSystems,numberOfSpecies),sysListRef(numberOfDisplacedSystems,numberOfSpecies)) + allocate(orbListI(numberOfDisplacedSystems),orbListII(numberOfDisplacedSystems)) + allocate(superMergedCoefficients(numberOfSpecies)) + allocate(superOverlapMatrix(numberOfSpecies), superMomentMatrix(numberOfSpecies,3),inverseOverlapMatrix(numberOfSpecies),molecularMomentMatrix(numberOfSpecies,3)) + allocate(franckCondonMatrix(numberOfSpecies),transitionDipoleMatrix(numberOfSpecies+1,3),refCurOverlapMatrix(numberOfSpecies),refCurMomentMatrix(numberOfSpecies,3)) + + arguments(1) = "NOCI-CONFIGURATIONCOEFFICIENTS" + ciCoefficients = Matrix_getFromFile(numberOfDisplacedSystems,numberOfDisplacedSystems,nociUnit,binary=.true.,arguments=arguments(1:1) ) + + arguments(1:1) = "NOCI-CONFIGURATIONENERGIES" + call Vector_getFromFile(numberOfDisplacedSystems, nociUnit, output=ciEnergies, binary=.true., arguments=arguments(1:1) ) + + arguments(1) = "MERGEDCOEFFICIENTS" + do speciesID=1, numberOfSpecies + numberOfContractions=molecularSystem_getTotalNumberOfContractions(speciesID) + dim2=max(MolecularSystem_getTotalNumberOfContractions(speciesID),MolecularSystem_getOcupationNumber(speciesID)) + arguments(2) = trim(MolecularSystem_instance%species(speciesID)%name) + auxCoefficients(speciesID) = Matrix_getFromFile(numberOfContractions,dim2,nociUnit,binary=.true.,arguments=arguments(1:2) ) + end do + + do sysI=1, numberOfDisplacedSystems + do speciesID=1, numberOfSpecies + numberOfContractions=molecularSystem_getTotalNumberOfContractions(speciesID) + write(auxString,*) sysI + arguments(1) = "SYSBASISLIST"//trim(auxString) + arguments(2) = trim(MolecularSystem_instance%species(speciesID)%name) + call Vector_getFromFileInteger(numberOfContractions, nociUnit, output=sysListRef(sysI,speciesID), binary=.true., arguments=arguments(1:2) ) + end do + end do + + close(nociUnit) + + !Create a super-mega molecular system + !Merge coefficients from NOCI calculation and reference system + + print *, "super-mega molecular system" + call MolecularSystem_mergeTwoSystems(superMergedMolecularSystem, this%mergedMolecularSystem, MolecularSystem_instance, & + orbListI(:),orbListII(:), reorder=.false.) + call MolecularSystem_showInformation(superMergedMolecularSystem) + call MolecularSystem_showParticlesInformation(superMergedMolecularSystem) + call MolecularSystem_showCartesianMatrix(superMergedMolecularSystem) + + call NonOrthogonalCI_mergeCoefficients(this%mergedCoefficients(:),auxCoefficients(:),& + this%mergedMolecularSystem,MolecularSystem_instance,superMergedMolecularSystem,& + orbListI(:),orbListII(:),superMergedCoefficients(:)) + + ! do speciesID=1, numberOfSpecies + ! print *, "superMergedCoefficients", speciesID + ! call Matrix_show(superMergedCoefficients(speciesID)) + ! end do + + !Fix basis list size + do speciesID=1, numberOfSpecies + ! print *, "orbListI", "speciesID", speciesID + ! call Vector_showInteger(orbListI(speciesID)) + do sysI=1, this%numberOfDisplacedSystems + call Vector_copyConstructorInteger(auxIVector,this%sysBasisList(sysI,speciesID)) + call Vector_constructorInteger(sysListCur(sysI,speciesID), MolecularSystem_getTotalNumberOfContractions(speciesID,superMergedMolecularSystem), 0) + do i=1, size(auxIVector%values) + if(orbListI(speciesID)%values(i) .eq. 0) cycle + sysListCur(sysI,speciesID)%values(i)=auxIVector%values(orbListI(speciesID)%values(i)) + end do + ! print *, "sysListCur", "sysI", sysI, "speciesID", speciesID + ! call Vector_showInteger(sysListCur(sysI,speciesID)) + end do + end do + + do speciesID=1, numberOfSpecies + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(1)) !not using the merged molecular systems + ! print *, "orbListII", "speciesID", speciesID + ! call Vector_showInteger(orbListII(speciesID)) + do sysII=1, numberOfDisplacedSystems + call Vector_copyConstructorInteger(auxIVector,sysListRef(sysII,speciesID)) + call Vector_constructorInteger(sysListRef(sysII,speciesID), MolecularSystem_getTotalNumberOfContractions(speciesID,superMergedMolecularSystem), 0) + do i=1, size(orbListII(speciesID)%values) + if(orbListII(speciesID)%values(i) .eq. 0) cycle + sysListRef(sysII,speciesID)%values(i)=auxIVector%values(orbListII(speciesID)%values(i)) + end do + ! print *, "sysListRef", "sysII", sysII, "speciesID", speciesID + ! call Vector_showInteger(sysListRef(sysII,speciesID)) + end do + end do + + ! if(CONTROL_instance%CI_STATES_TO_PRINT .eq. 0) return + + ! numberOfSpecies=molecularSystem_instance%numberOfQuantumSpecies + + + print *, "" + print *, "Computing overlap and moment integrals for the super-mega system..." + print *, "" + do speciesID = 1, numberOfSpecies + call DirectIntegralManager_getOverlapIntegrals(superMergedMolecularSystem,speciesID,superOverlapMatrix(speciesID)) + call DirectIntegralManager_getMomentIntegrals(superMergedMolecularSystem,speciesID,1,superMomentMatrix(speciesID,1)) + call DirectIntegralManager_getMomentIntegrals(superMergedMolecularSystem,speciesID,2,superMomentMatrix(speciesID,2)) + call DirectIntegralManager_getMomentIntegrals(superMergedMolecularSystem,speciesID,3,superMomentMatrix(speciesID,3)) + end do + !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for supermolecular 1-body integrals : ", omp_get_wtime() - timeA ," (s)" + !$ timeA = omp_get_wtime() + + print *, "" + print *, "Self overlap matrices for the supermegaposed systems..." + print *, "" + + do speciesID=1, numberOfSpecies + call Matrix_constructor(refCurOverlapMatrix(speciesID), int(this%numberOfDisplacedSystems,8), & + int(numberOfDisplacedSystems,8), 1.0_8) + end do + !!Fill the merged density matrix + !!"Non Diagonal" terms - system pairs + do sysI=1, numberOfDisplacedSystems !computed + do sysII=1, numberOfDisplacedSystems !reference + ! if( abs(this%configurationOverlapMatrix%values(sysI,sysII)) .lt. CONTROL_instance%CONFIGURATION_OVERLAP_THRESHOLD ) cycle + !!Compute molecular overlap matrix and its inverse + do speciesID=1, numberOfSpecies + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) !not using the merged molecular systems + particlesPerOrbital=MolecularSystem_getEta(speciesID,this%molecularSystems(sysI)) + call Matrix_constructor(molecularOverlapMatrix, int(occupationNumber,8), int(occupationNumber,8), 0.0_8 ) + ! call Matrix_constructor(inverseOverlapMatrix, int(occupationNumber,8), int(occupationNumber,8), 0.0_8 ) + + do mu=1, MolecularSystem_getTotalNumberOfContractions(speciesID,superMergedMolecularSystem) !sysI + if(sysListRef(sysI,speciesID)%values(mu) .eq. 0) cycle + do nu=1, MolecularSystem_getTotalNumberOfContractions(speciesID,superMergedMolecularSystem) !sysII + if(sysListRef(sysII,speciesID)%values(nu) .eq. 0) cycle + do i = 1 , occupationNumber + ii=occupationNumber*(sysI-1)+i+MolecularSystem_getOcupationNumber(speciesID,superMergedMolecularSystem)/2 + do j = 1 , occupationNumber + jj=occupationNumber*(sysII-1)+j+MolecularSystem_getOcupationNumber(speciesID,superMergedMolecularSystem)/2 + ! print *, "i, j, mu, nu, coefI, coefII, overlap", i,j,mu,nu,superMergedCoefficients(speciesID)%values(mu,ii),& + ! superMergedCoefficients(speciesID)%values(nu,jj),& + ! superOverlapMatrix(speciesID)%values(mu,nu) + molecularOverlapMatrix%values(i,j)=molecularOverlapMatrix%values(i,j)+& + superMergedCoefficients(speciesID)%values(mu,ii)*& + superMergedCoefficients(speciesID)%values(nu,jj)*& + superOverlapMatrix(speciesID)%values(mu,nu) + end do + end do + end do + end do + if(occupationNumber .ne. 0) then + ! inverseOverlapMatrix=Matrix_inverse(molecularOverlapMatrix) + ! print *, "inverseOverlapMatrices sysI, sysII", speciesID, sysI, sysII + ! call Matrix_show(inverseOverlapMatrices(speciesID)) + call Matrix_getDeterminant(molecularOverlapMatrix,overlapDeterminant,method="LU") + ! print *, "OverlapDeterminantLU speciesID, sysI, sysII", speciesID, sysI, sysII, overlapDeterminant + else + overlapDeterminant=1.0 + end if + refCurOverlapMatrix(speciesID)%values(sysI,sysII)=refCurOverlapMatrix(speciesID)%values(sysI,sysII)*overlapDeterminant**particlesPerOrbital + end do + + end do + end do + + do speciesID=1, numberOfSpecies + print *, "Reference Overlap Matrix for", speciesID + call Matrix_show(refCurOverlapMatrix(speciesID)) + end do + + print *, "" + print *, "Building Franck-Condon matrices for the superposed systems..." + print *, "" + + do speciesID=1, numberOfSpecies + call Matrix_constructor(refCurOverlapMatrix(speciesID), int(this%numberOfDisplacedSystems,8), & + int(numberOfDisplacedSystems,8), 1.0_8) + do k=1,3 + call Matrix_constructor(refCurMomentMatrix(speciesID,k), int(this%numberOfDisplacedSystems,8), & + int(numberOfDisplacedSystems,8), 0.0_8) + end do + end do + call Matrix_constructor(refCurTotalOverlap, int(this%numberOfDisplacedSystems,8), & + int(numberOfDisplacedSystems,8), 1.0_8) + + !!Fill the merged density matrix + !!"Non Diagonal" terms - system pairs + do sysI=1, this%numberOfDisplacedSystems !computed + do sysII=1, numberOfDisplacedSystems !reference + ! if( abs(this%configurationOverlapMatrix%values(sysI,sysII)) .lt. CONTROL_instance%CONFIGURATION_OVERLAP_THRESHOLD ) cycle + !!Compute molecular overlap matrix and its inverse + do speciesID=1, numberOfSpecies + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) !not using the merged molecular systems + particlesPerOrbital=MolecularSystem_getEta(speciesID,this%molecularSystems(sysI)) + call Matrix_constructor(molecularOverlapMatrix, int(occupationNumber,8), int(occupationNumber,8), 0.0_8 ) + do k=1,3 + call Matrix_constructor(molecularMomentMatrix(speciesID,k), int(occupationNumber,8), int(occupationNumber,8), 0.0_8 ) + end do + + do mu=1, MolecularSystem_getTotalNumberOfContractions(speciesID,superMergedMolecularSystem) !sysI + if(sysListCur(sysI,speciesID)%values(mu) .eq. 0) cycle + do nu=1, MolecularSystem_getTotalNumberOfContractions(speciesID,superMergedMolecularSystem) !sysII + if(sysListRef(sysII,speciesID)%values(nu) .eq. 0) cycle + do i = 1 , occupationNumber + ii=occupationNumber*(sysI-1)+i + do j = 1 , occupationNumber + jj=occupationNumber*(sysII-1)+j+MolecularSystem_getOcupationNumber(speciesID,superMergedMolecularSystem)/2 + ! print *, "i, j, mu, nu, coefI, coefII, overlap", i,j,mu,nu,superMergedCoefficients(speciesID)%values(mu,ii),& + ! superMergedCoefficients(speciesID)%values(nu,jj),& + ! superOverlapMatrix(speciesID)%values(mu,nu) + molecularOverlapMatrix%values(i,j)=molecularOverlapMatrix%values(i,j)+& + superMergedCoefficients(speciesID)%values(mu,ii)*& + superMergedCoefficients(speciesID)%values(nu,jj)*& + superOverlapMatrix(speciesID)%values(mu,nu) + do k=1,3 + molecularMomentMatrix(speciesID,k)%values(i,j)=molecularMomentMatrix(speciesID,k)%values(i,j)+& + superMergedCoefficients(speciesID)%values(mu,ii)*& + superMergedCoefficients(speciesID)%values(nu,jj)*& + superMomentMatrix(speciesID,k)%values(mu,nu) + end do + end do + end do + end do + end do + if(occupationNumber .ne. 0) then + inverseOverlapMatrix(speciesID)=Matrix_inverse(molecularOverlapMatrix) + ! print *, "inverseOverlapMatrices sysI, sysII", speciesID, sysI, sysII + ! call Matrix_show(inverseOverlapMatrices(speciesID)) + call Matrix_getDeterminant(molecularOverlapMatrix,overlapDeterminant,method="LU") + ! print *, "OverlapDeterminantLU speciesID, sysI, sysII", speciesID, sysI, sysII, overlapDeterminant + refCurOverlapMatrix(speciesID)%values(sysI,sysII)=overlapDeterminant**particlesPerOrbital + else + overlapDeterminant=1.0 + end if + refCurTotalOverlap%values(sysI,sysII)=refCurTotalOverlap%values(sysI,sysII)*refCurOverlapMatrix(speciesID)%values(sysI,sysII) + end do + + do speciesID=1, numberOfSpecies + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(sysI)) !not using the merged molecular systems + particlesPerOrbital=MolecularSystem_getEta(speciesID,this%molecularSystems(sysI)) + do i = 1 , occupationNumber + do j = 1 , occupationNumber + do k=1,3 + refCurMomentMatrix(speciesID,k)%values(sysI,sysII)=refCurMomentMatrix(speciesID,k)%values(sysI,sysII)+& + molecularMomentMatrix(speciesID,k)%values(i,j)*& + inverseOverlapMatrix(speciesID)%values(j,i) + end do + end do + end do + do k=1,3 + refCurMomentMatrix(speciesID,k)%values(sysI,sysII)=refCurMomentMatrix(speciesID,k)%values(sysI,sysII)*refCurTotalOverlap%values(sysI,sysII)*particlesPerOrbital + end do + end do + end do + end do + + do speciesID=1, numberOfSpecies + print *, "refCurOverlapMatrix(speciesID)", speciesID + call Matrix_show(refCurOverlapMatrix(speciesID)) + call Matrix_constructor(franckCondonMatrix(speciesID), int(CONTROL_instance%CI_STATES_TO_PRINT,8), int(CONTROL_instance%CI_STATES_TO_PRINT,8), 0.0_8) + end do + + !+1 For point charges + do speciesID=1, numberOfSpecies+1 + do k=1,3 + call Matrix_constructor(transitionDipoleMatrix(speciesID,k), int(CONTROL_instance%CI_STATES_TO_PRINT,8), int(CONTROL_instance%CI_STATES_TO_PRINT,8), 0.0_8) + end do + end do + + do stateII=1, CONTROL_instance%CI_STATES_TO_PRINT + print *, "Reference state:", stateII + do stateI=1, CONTROL_instance%CI_STATES_TO_PRINT + print *, " current state:", stateI + do speciesID=1, numberOfSpecies + occupationNumber=MolecularSystem_getOcupationNumber(speciesID,this%molecularSystems(1)) !not using the merged molecular systems + print *, "occupationNumber", occupationNumber + particlesPerOrbital=MolecularSystem_getEta(speciesID,this%molecularSystems(1)) + trololo=0 + do sysI=1, this%numberOfDisplacedSystems !computed + do sysII=1, numberOfDisplacedSystems !reference + do i = 1 , occupationNumber + do j = 1 , occupationNumber + trololo = trololo + & + inverseOverlapMatrix(speciesID)%values(j,i)*& + this%configurationCoefficients%values(sysI,stateI)*& + ciCoefficients%values(sysII,stateII)*& !!reference + refCurOverlapMatrix(speciesID)%values(sysI,sysII)*& + particlesPerOrbital + end do + end do + ! refCurTotalOverlap%values(sysI,sysII)*& + ! franckCondonMatrix(speciesID)%values(stateI,stateII)+& + + do k=1,3 + transitionDipoleMatrix(speciesID,k)%values(stateI,stateII) = transitionDipoleMatrix(speciesID,k)%values(stateI,stateII) + & + molecularsystem_getcharge( speciesID )*& + this%configurationCoefficients%values(sysI,stateI)*& + ciCoefficients%values(sysII,stateII)*& !!reference + refCurMomentMatrix(speciesID,k)%values(sysI,sysII) + end do + + end do + end do + print *, "speciesID", speciesID, "trololo", trololo + franckCondonMatrix(speciesID)%values(stateI,stateII)=trololo + franckCondonMatrix(speciesID)%values(stateI,stateII)=franckCondonMatrix(speciesID)%values(stateI,stateII)/(occupationNumber*particlesPerOrbital) + print *, " F.C. factor for ", molecularSystem_getNameOfSpecies(speciesID),& + franckCondonMatrix(speciesID)%values(stateI,stateII) + end do + do sysI=1, this%numberOfDisplacedSystems !computed + do sysII=1, numberOfDisplacedSystems !reference + do k=1,3 + transitionDipoleMatrix(numberOfSpecies+1,k)%values(stateI,stateII) = transitionDipoleMatrix(numberOfSpecies+1,k)%values(stateI,stateII) + & + pointchargesdipole(k)*& + this%configurationCoefficients%values(sysI,stateI)*& + ciCoefficients%values(sysII,stateII)*& !!reference + refCurTotalOverlap%values(sysI,sysII) + end do + end do + end do + ! trololo=1 + ! do speciesID=1, numberOfSpecies + ! trololo=trololo*franckCondonMatrix(speciesID)%values(stateI,stateII) + ! end do + ! print *, " F.C. factor product ", trololo + ! trololo=0 + ! do speciesID=1, numberOfSpecies + ! trololo=trololo+franckCondonMatrix(speciesID)%values(stateI,stateII) + ! end do + ! print *, " F.C. factor sum ", trololo + ! trololo=0 + ! do sysI=1, this%numberOfDisplacedSystems !computed + ! do sysII=1, numberOfDisplacedSystems !reference + ! trololo = trololo + & + ! this%configurationCoefficients%values(sysI,stateI)*& + ! ciCoefficients%values(sysII,stateII)*& !!reference + ! refCurTotalOverlap%values(sysI,sysII) + ! end do + ! end do + ! print *, " total overlap ", trololo + end do + end do + + print *, "Dipole approximation spectrum" + do stateII=1, CONTROL_instance%CI_STATES_TO_PRINT + print *, "Reference state:", stateII + do stateI=1, CONTROL_instance%CI_STATES_TO_PRINT + trolololo=0 + print *, "current state:", stateI + do speciesID=1, numberOfSpecies + do k=1,3 + trolololo(k)=trolololo(k)+transitionDipoleMatrix(speciesID,k)%values(stateI,stateII) + end do + print *, " T.D. integrals for ", molecularSystem_getNameOfSpecies(speciesID),& + transitionDipoleMatrix(speciesID,1)%values(stateI,stateII),& + transitionDipoleMatrix(speciesID,2)%values(stateI,stateII),& + transitionDipoleMatrix(speciesID,3)%values(stateI,stateII) + end do + do k=1,3 + trolololo(k)=trolololo(k)+transitionDipoleMatrix(numberOfSpecies+1,k)%values(stateI,stateII) + end do + print *, " T.D. integrals point charges ", & + transitionDipoleMatrix(numberOfSpecies+1,1)%values(stateI,stateII),& + transitionDipoleMatrix(numberOfSpecies+1,2)%values(stateI,stateII),& + transitionDipoleMatrix(numberOfSpecies+1,3)%values(stateI,stateII) + print *, "energy dif", ciEnergies%values(stateII)-this%statesEigenvalues%values(stateI), "total components", trolololo(1:3) ,"intensity", sqrt(sum(trolololo(1:3)**2)) + end do + end do + + close(densUnit) + + deallocate(auxCoefficients,& + sysListCur,sysListRef,& + orbListI,orbListII,& + superMergedCoefficients,& + superOverlapMatrix,& + franckCondonMatrix) + + end subroutine NonOrthogonalCI_computeFranckCondon + end module NonOrthogonalCI_ diff --git a/src/PT/PropagatorTheory.f90 b/src/PT/PropagatorTheory.f90 index c5b882da..f7616aa8 100644 --- a/src/PT/PropagatorTheory.f90 +++ b/src/PT/PropagatorTheory.f90 @@ -274,11 +274,11 @@ subroutine PropagatorTheory_show() write (6,"(T10,A40)") "-------------------------------------" end if - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE") then - species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE(1) ) - do z = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if (CONTROL_instance%IONIZE_SPECIE(z) /= "NONE" ) then - species2ID= MolecularSystem_getSpecieID(CONTROL_instance%IONIZE_SPECIE(z)) + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE") then + species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES(1) ) + do z = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if (CONTROL_instance%IONIZE_SPECIES(z) /= "NONE" ) then + species2ID= MolecularSystem_getSpecieID(CONTROL_instance%IONIZE_SPECIES(z)) end if end do else @@ -385,11 +385,11 @@ subroutine PropagatorTheory_show() write (*,"(T10,A60)") "LOWDIN implementation: J. Chem. Phys. 138, 194108 (2013)" write (*,"(T10,A64)") "---------------------------------------------------------------" - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE") then - species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE(1) ) - do z = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if (CONTROL_instance%IONIZE_SPECIE(z) /= "NONE" ) then - species2ID= MolecularSystem_getSpecieID(CONTROL_instance%IONIZE_SPECIE(z)) + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE") then + species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES(1) ) + do z = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if (CONTROL_instance%IONIZE_SPECIES(z) /= "NONE" ) then + species2ID= MolecularSystem_getSpecieID(CONTROL_instance%IONIZE_SPECIES(z)) end if end do else @@ -521,8 +521,8 @@ subroutine PropagatorTheory_secondOrderCorrection() !!! Defining for which species the correction will be applied -! if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE") then -! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE(1) ) +! if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE") then +! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES(1) ) ! species2ID= species1ID ! m=1 ! else @@ -531,13 +531,13 @@ subroutine PropagatorTheory_secondOrderCorrection() ! m = species2ID ! end if - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE") then - species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE(1) ) + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE") then + species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES(1) ) zz = 0 - do z = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if (CONTROL_instance%IONIZE_SPECIE(z) /= "NONE" ) then + do z = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if (CONTROL_instance%IONIZE_SPECIES(z) /= "NONE" ) then zz = zz + 1 - species2ID= MolecularSystem_getSpecieID(CONTROL_instance%IONIZE_SPECIE(z)) + species2ID= MolecularSystem_getSpecieID(CONTROL_instance%IONIZE_SPECIES(z)) end if end do m = zz @@ -604,16 +604,16 @@ subroutine PropagatorTheory_secondOrderCorrection() !!! Defining the number of orbitals !!! Insert a parameter for the else option if (CONTROL_instance%PT_TRANSITION_OPERATOR.or.CONTROL_instance%PT_JUST_ONE_ORBITAL) then - PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO - PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO + PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) + PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) n = 1 - !else if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE".and.CONTROL_instance%IONIZE_MO /= 0) then + !else if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE".and.CONTROL_instance%IONIZE_MO(1) /= 0) then !PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 - !PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO + !PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO(1) !n = PropagatorTheory_instance%virtualBoundary-PropagatorTheory_instance%occupationBoundary+1 - else if ( CONTROL_instance%IONIZE_MO /= 0) then - PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO - PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO + else if ( CONTROL_instance%IONIZE_MO(1) /= 0) then + PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) + PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) n = 1 else PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 @@ -641,7 +641,7 @@ subroutine PropagatorTheory_secondOrderCorrection() if (CONTROL_instance%PT_TRANSITION_OPERATOR) then - occupationsOfSpeciesA%values(CONTROL_instance%IONIZE_MO)=CONTROL_instance%MO_FRACTION_OCCUPATION + occupationsOfSpeciesA%values(CONTROL_instance%IONIZE_MO(1))=CONTROL_instance%MO_FRACTION_OCCUPATION(1) end if @@ -1347,8 +1347,8 @@ subroutine PropagatorTheory_thirdOrderCorrection5() ! Defining for which species the correction will be applied - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE") then - species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE(1) ) + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE") then + species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES(1) ) species2ID= species1ID m=1 else @@ -1563,12 +1563,12 @@ subroutine PropagatorTheory_thirdOrderCorrection5() ! Defining the number of orbitals !!! Insert a parameter for the else option if (CONTROL_instance%PT_JUST_ONE_ORBITAL) then - PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO - PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO + PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) + PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) n = 1 - else if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE".and.CONTROL_instance%IONIZE_MO /= 0) then + else if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE".and.CONTROL_instance%IONIZE_MO(1) /= 0) then PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 - PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO + PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO(1) n = PropagatorTheory_instance%virtualBoundary-PropagatorTheory_instance%occupationBoundary+1 else PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 @@ -4861,7 +4861,7 @@ end module PropagatorTheory_ ! end if ! if (CONTROL_instance%PT_TRANSITION_OPERATOR.or.CONTROL_instance%PT_JUST_ONE_ORBITAL) then - ! specie1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) + ! specie1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! specie2ID= specie1ID ! else ! specie1ID=1 @@ -4880,8 +4880,8 @@ end module PropagatorTheory_ ! lambda = MolecularSystem_getLambda( i ) ! if (CONTROL_instance%PT_TRANSITION_OPERATOR.or.CONTROL_instance%PT_JUST_ONE_ORBITAL) then - ! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO - ! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO + ! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) + ! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! else ! PropagatorTheory_instance%virtualBoundary= occupationNumber ! PropagatorTheory_instance%occupationBoundary = 1 @@ -4905,8 +4905,8 @@ end module PropagatorTheory_ ! if (CONTROL_instance%PT_TRANSITION_OPERATOR) then ! call Vector_constructor(occupations,occupationNumber) ! do k=1, occupationNumber - ! if (k==CONTROL_instance%IONIZE_MO) then - ! occupations%values(k)=CONTROL_instance%MO_FRACTION_OCCUPATION + ! if (k==CONTROL_instance%IONIZE_MO(1)) then + ! occupations%values(k)=CONTROL_instance%MO_FRACTION_OCCUPATION(1) ! else ! occupations%values(k)=1.0_8 ! end if @@ -4998,7 +4998,7 @@ end module PropagatorTheory_ ! end do ! end do - ! if (CONTROL_instance%PT_TRANSITION_OPERATOR .and. nameOfSpecie == CONTROL_instance%IONIZE_SPECIE .and. p==CONTROL_instance%IONIZE_MO) then + ! if (CONTROL_instance%PT_TRANSITION_OPERATOR .and. nameOfSpecie == CONTROL_instance%IONIZE_SPECIES .and. p==CONTROL_instance%IONIZE_MO(1)) then ! do a=1, occupationNumber ! do s=occupationNumber+1, numberOfContractions ! id=id+1 @@ -5046,8 +5046,8 @@ end module PropagatorTheory_ ! if (nameOfOtherSpecie == CONTROL_instance%IONIZE_SPECIE) then ! call Vector_constructor(occupationsOfOtherSpecie,occupationNumberOfOtherSpecie) ! do k=1, occupationNumberOfOtherSpecie - ! if (k==CONTROL_instance%IONIZE_MO) then - ! occupationsOfOtherSpecie%values(k)=CONTROL_instance%MO_FRACTION_OCCUPATION + ! if (k==CONTROL_instance%IONIZE_MO(1)) then + ! occupationsOfOtherSpecie%values(k)=CONTROL_instance%MO_FRACTION_OCCUPATION(1) ! else ! occupationsOfOtherSpecie%values(k)=1.0_8 ! end if @@ -5135,7 +5135,7 @@ end module PropagatorTheory_ ! end do ! end do - ! if (CONTROL_instance%PT_TRANSITION_OPERATOR .and. nameOfSpecie == CONTROL_instance%IONIZE_SPECIE .and. p==CONTROL_instance%IONIZE_MO) then + ! if (CONTROL_instance%PT_TRANSITION_OPERATOR .and. nameOfSpecie == CONTROL_instance%IONIZE_SPECIES .and. p==CONTROL_instance%IONIZE_MO(1)) then ! do a=1, occupationNumberOfOtherSpecie ! do s=occupationNumberOfOtherSpecie+1, numberOfContractionsOfOtherSpecie @@ -5298,7 +5298,7 @@ end module PropagatorTheory_ ! print *,"******************************************************************" ! print *,"BEGINNING OF SECOND ORDER ELECTRON-NUCLEAR PROPAGATOR CALCULATIONS" - ! speciesID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) + ! speciesID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! nameOfSpecies= trim( MolecularSystem_getNameOfSpecie( speciesID ) ) @@ -5315,8 +5315,8 @@ end module PropagatorTheory_ ! lambdaOfSpecies = MolecularSystem_getLambda( speciesID ) ! if ( CONTROL_instance%PT_JUST_ONE_ORBITAL) then - ! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO - ! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO + ! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) + ! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! else ! PropagatorTheory_instance%virtualBoundary= occupationNumberOfSpecies ! PropagatorTheory_instance%occupationBoundary = 1 @@ -5579,7 +5579,7 @@ end module PropagatorTheory_ ! do pa = 1, HamiltonianSize - ! if ( pa==CONTROL_instance%IONIZE_MO ) then + ! if ( pa==CONTROL_instance%IONIZE_MO(1) ) then ! lastEigenvector%values(pa) = 1.0_8 @@ -5591,7 +5591,7 @@ end module PropagatorTheory_ ! end do - ! lastEigenvalue = eigenValuesOfSpecies%values(CONTROL_instance%IONIZE_MO) + ! lastEigenvalue = eigenValuesOfSpecies%values(CONTROL_instance%IONIZE_MO(1)) ! print *,lastEigenvector%values @@ -5715,7 +5715,7 @@ end module PropagatorTheory_ ! ! Storing of corrections - ! PropagatorTheory_instance%energyCorrectionsOfSecondOrder%values(speciesID,4*CONTROL_instance%IONIZE_MO)= & + ! PropagatorTheory_instance%energyCorrectionsOfSecondOrder%values(speciesID,4*CONTROL_instance%IONIZE_MO(1))= & ! 27.211396_8 * lastEigenvalue ! !!! DAVIDSON ALGORYTHM ENDS @@ -5800,7 +5800,7 @@ end module PropagatorTheory_ ! print *,"******************************************************************" ! print *,"BEGINNING OF SECOND ORDER ELECTRON-NUCLEAR PROPAGATOR CALCULATIONS" -! speciesID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) +! speciesID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! nameOfSpecies= trim( MolecularSystem_getNameOfSpecie( speciesID ) ) ! chargeOfSpecies = MolecularSystem_getCharge( speciesID ) ! eigenValuesOfSpecies = MolecularSystem_getEigenValues( speciesID ) @@ -6563,8 +6563,8 @@ end module PropagatorTheory_ ! print *,"idaHf:",idaHf,"idfHf:",idfHf,"HamiltonianSize:",HamiltonianSize ! ! if ( CONTROL_instance%PT_JUST_ONE_ORBITAL) then -! ! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO -! ! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO +! ! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) +! ! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! ! else ! ! PropagatorTheory_instance%virtualBoundary= occupationNumberOfSpecies ! ! PropagatorTheory_instance%occupationBoundary = 1 @@ -6581,9 +6581,9 @@ end module PropagatorTheory_ ! lastEigenvector%values = 0.0_8 -! lastEigenvector%values(CONTROL_instance%IONIZE_MO) = 1.0_8 +! lastEigenvector%values(CONTROL_instance%IONIZE_MO(1)) = 1.0_8 -! lastEigenvalue = eigenValuesOfSpecies%values(CONTROL_instance%IONIZE_MO) +! lastEigenvalue = eigenValuesOfSpecies%values(CONTROL_instance%IONIZE_MO(1)) ! print *,"Initial lastEigenvector values:" ! call Vector_show(lastEigenvector) @@ -6595,10 +6595,10 @@ end module PropagatorTheory_ ! call Matrix_eigen(superHamiltonian, superEigenvalues, superEigenvectors, SYMMETRIC) -! print *,"koopmans:",eigenvaluesOfSpecies%values(CONTROL_instance%IONIZE_MO) +! print *,"koopmans:",eigenvaluesOfSpecies%values(CONTROL_instance%IONIZE_MO(1)) ! ! call Vector_show(superEigenvalues) -! ! print *,"superEigenvalue:",superEigenvalues%values(CONTROL_instance%IONIZE_MO) +! ! print *,"superEigenvalue:",superEigenvalues%values(CONTROL_instance%IONIZE_MO(1)) ! do k = 1, HamiltonianSize ! dot=dot_product(superEigenvectors%values(:,k),lastEigenvector%values) @@ -6862,7 +6862,7 @@ end module PropagatorTheory_ ! ! Storing of corrections -! PropagatorTheory_instance%energyCorrectionsOfSecondOrder%values(speciesID,4*CONTROL_instance%IONIZE_MO)= & +! PropagatorTheory_instance%energyCorrectionsOfSecondOrder%values(speciesID,4*CONTROL_instance%IONIZE_MO(1))= & ! 27.211396_8 * lastEigenvalue ! ! PropagatorTheory_instance%energyCorrectionsOfSecondOrder%values(i,(4*m-1))= & @@ -6951,8 +6951,8 @@ end module PropagatorTheory_ ! ! !!! Defining for which species the correction will be applied ! -! if (CONTROL_instance%IONIZE_SPECIE /= "NONE") then -! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) +! if (CONTROL_instance%IONIZE_SPECIES /= "NONE") then +! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! species2ID= species1ID ! m=1 ! else @@ -6987,12 +6987,12 @@ end module PropagatorTheory_ ! !!! Defining the number of orbitals !!! Insert a parameter for the else option ! ! if (CONTROL_instance%PT_JUST_ONE_ORBITAL) then -! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO -! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) +! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! n = 1 -! else if (CONTROL_instance%IONIZE_SPECIE /= "NONE".and.CONTROL_instance%IONIZE_MO /= 0) then +! else if (CONTROL_instance%IONIZE_SPECIES /= "NONE".and.CONTROL_instance%IONIZE_MO(1) /= 0) then ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 -! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO(1) ! n = PropagatorTheory_instance%virtualBoundary-PropagatorTheory_instance%occupationBoundary+1 ! else ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 @@ -8394,8 +8394,8 @@ end module PropagatorTheory_ ! ! !!! Defining for which species the correction will be applied ! -! if (CONTROL_instance%IONIZE_SPECIE /= "NONE") then -! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) +! if (CONTROL_instance%IONIZE_SPECIES /= "NONE") then +! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! species2ID= species1ID ! m=1 ! else @@ -8430,12 +8430,12 @@ end module PropagatorTheory_ ! ! Defining the number of orbitals !!! Insert a parameter for the else option ! ! if (CONTROL_instance%PT_JUST_ONE_ORBITAL) then -! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO -! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) +! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! n = 1 -! else if (CONTROL_instance%IONIZE_SPECIE /= "NONE".and.CONTROL_instance%IONIZE_MO /= 0) then +! else if (CONTROL_instance%IONIZE_SPECIES /= "NONE".and.CONTROL_instance%IONIZE_MO(1) /= 0) then ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 -! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO(1) ! n = PropagatorTheory_instance%virtualBoundary-PropagatorTheory_instance%occupationBoundary+1 ! else ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 @@ -10546,8 +10546,8 @@ end module PropagatorTheory_ ! ! !!! Defining for which species the correction will be applied ! -! if (CONTROL_instance%IONIZE_SPECIE /= "NONE") then -! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) +! if (CONTROL_instance%IONIZE_SPECIES /= "NONE") then +! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! species2ID= species1ID ! m=1 ! else @@ -10617,12 +10617,12 @@ end module PropagatorTheory_ ! ! Defining the number of orbitals !!! Insert a parameter for the else option ! ! if (CONTROL_instance%PT_JUST_ONE_ORBITAL) then -! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO -! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) +! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! n = 1 -! else if (CONTROL_instance%IONIZE_SPECIE /= "NONE".and.CONTROL_instance%IONIZE_MO /= 0) then +! else if (CONTROL_instance%IONIZE_SPECIES /= "NONE".and.CONTROL_instance%IONIZE_MO(1) /= 0) then ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 -! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO(1) ! n = PropagatorTheory_instance%virtualBoundary-PropagatorTheory_instance%occupationBoundary+1 ! else ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 @@ -12723,8 +12723,8 @@ end module PropagatorTheory_ ! ! ! Defining for which species the correction will be applied ! -! if (CONTROL_instance%IONIZE_SPECIE /= "NONE") then -! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIE ) +! if (CONTROL_instance%IONIZE_SPECIES /= "NONE") then +! species1ID = MolecularSystem_getSpecieID( nameOfSpecie=CONTROL_instance%IONIZE_SPECIES ) ! species2ID= species1ID ! m=1 ! else @@ -12807,12 +12807,12 @@ end module PropagatorTheory_ ! ! Defining the number of orbitals !!! Insert a parameter for the else option ! ! if (CONTROL_instance%PT_JUST_ONE_ORBITAL) then -! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO -! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%virtualBoundary=CONTROL_instance%IONIZE_MO(1) +! PropagatorTheory_instance%occupationBoundary=CONTROL_instance%IONIZE_MO(1) ! n = 1 -! else if (CONTROL_instance%IONIZE_SPECIE /= "NONE".and.CONTROL_instance%IONIZE_MO /= 0) then +! else if (CONTROL_instance%IONIZE_SPECIES /= "NONE".and.CONTROL_instance%IONIZE_MO(1) /= 0) then ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 -! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO +! PropagatorTheory_instance%occupationBoundary = CONTROL_instance%IONIZE_MO(1) ! n = PropagatorTheory_instance%virtualBoundary-PropagatorTheory_instance%occupationBoundary+1 ! else ! PropagatorTheory_instance%virtualBoundary = occupationNumberOfSpeciesA + 1 diff --git a/src/core/CONTROL.f90 b/src/core/CONTROL.f90 index 946459bd..ad29e403 100644 --- a/src/core/CONTROL.f90 +++ b/src/core/CONTROL.f90 @@ -57,6 +57,7 @@ module CONTROL_ real(8) :: DIIS_SWITCH_THRESHOLD_BKP real(8) :: ELECTRONIC_LEVEL_SHIFTING real(8) :: NONELECTRONIC_LEVEL_SHIFTING + real(8) :: EXCHANGE_ORBITAL_THRESHOLD real(8) :: WAVE_FUNCTION_SCALE integer :: SCF_NONELECTRONIC_MAX_ITERATIONS integer :: SCF_ELECTRONIC_MAX_ITERATIONS @@ -332,10 +333,10 @@ module CONTROL_ !!***************************************************** !! Miscelaneous Options !! - real(8) :: MO_FRACTION_OCCUPATION - integer :: IONIZE_MO - character(50) :: IONIZE_SPECIE(10) - character(50) :: EXCITE_SPECIE + real(8) :: MO_FRACTION_OCCUPATION(10) + integer :: IONIZE_MO(10) + character(50) :: IONIZE_SPECIES(10) + character(50) :: EXCITE_SPECIES integer :: NUMBER_OF_CORES !!***************************************************** @@ -395,6 +396,7 @@ module CONTROL_ real(8) :: LowdinParameters_diisSwitchThreshold_bkp real(8) :: LowdinParameters_electronicLevelShifting real(8) :: LowdinParameters_nonelectronicLevelShifting + real(8) :: LowdinParameters_exchangeOrbitalThreshold real(8) :: LowdinParameters_waveFunctionScale integer :: LowdinParameters_scfNonelectronicMaxIterations integer :: LowdinParameters_scfElectronicMaxIterations @@ -669,10 +671,10 @@ module CONTROL_ !!***************************************************** !! Miscelaneous Options !! - real(8) :: LowdinParameters_MOFractionOccupation - integer :: LowdinParameters_ionizeMO - character(50) :: LowdinParameters_ionizeSpecie(10) - character(50) :: LowdinParameters_exciteSpecie + real(8) :: LowdinParameters_MOFractionOccupation(10) + integer :: LowdinParameters_ionizeMO(10) + character(50) :: LowdinParameters_ionizeSpecies(10) + character(50) :: LowdinParameters_exciteSpecies !!***************************************************** !! Integrals transformation options @@ -731,6 +733,7 @@ module CONTROL_ LowdinParameters_diisSwitchThreshold_bkp,& LowdinParameters_electronicLevelShifting,& LowdinParameters_nonelectronicLevelShifting,& + LowdinParameters_exchangeOrbitalThreshold,& LowdinParameters_waveFunctionScale,& LowdinParameters_scfNonelectronicMaxIterations,& LowdinParameters_scfElectronicMaxIterations,& @@ -1006,8 +1009,8 @@ module CONTROL_ !! LowdinParameters_MOFractionOccupation,& LowdinParameters_ionizeMO,& - LowdinParameters_ionizeSpecie,& - LowdinParameters_exciteSpecie,& + LowdinParameters_ionizeSpecies,& + LowdinParameters_exciteSpecies,& !!***************************************************** !! Integrals transformation options @@ -1091,6 +1094,7 @@ subroutine CONTROL_start() LowdinParameters_diisSwitchThreshold_bkp = 0.5 LowdinParameters_electronicLevelShifting = 0.0 LowdinParameters_nonelectronicLevelShifting = 0.0 + LowdinParameters_exchangeOrbitalThreshold = 0.8 LowdinParameters_waveFunctionScale = 1000.0 LowdinParameters_scfNonelectronicMaxIterations = 50 LowdinParameters_scfElectronicMaxIterations = 50 @@ -1366,9 +1370,8 @@ subroutine CONTROL_start() !! LowdinParameters_MOFractionOccupation = 1.0_8 LowdinParameters_ionizeMO = 0 - LowdinParameters_ionizeSpecie = "NONE" - LowdinParameters_exciteSpecie = "NONE" - LowdinParameters_exciteSpecie = "NONE" + LowdinParameters_ionizeSpecies = "NONE" + LowdinParameters_exciteSpecies = "NONE" !!***************************************************** !! Integrals transformation options @@ -1430,6 +1433,7 @@ subroutine CONTROL_start() CONTROL_instance%DIIS_SWITCH_THRESHOLD_BKP = 0.5 CONTROL_instance%ELECTRONIC_LEVEL_SHIFTING = 0.0 CONTROL_instance%NONELECTRONIC_LEVEL_SHIFTING = 0.0 + CONTROL_instance%EXCHANGE_ORBITAL_THRESHOLD = 0.8 CONTROL_instance%WAVE_FUNCTION_SCALE = 1000.0 CONTROL_instance%SCF_NONELECTRONIC_MAX_ITERATIONS = 50 CONTROL_instance%SCF_ELECTRONIC_MAX_ITERATIONS = 50 @@ -1702,8 +1706,8 @@ subroutine CONTROL_start() !! CONTROL_instance%MO_FRACTION_OCCUPATION = 1.0_8 CONTROL_instance%IONIZE_MO = 0 - CONTROL_instance%IONIZE_SPECIE = "NONE" - CONTROL_instance%EXCITE_SPECIE = "NONE" + CONTROL_instance%IONIZE_SPECIES = "NONE" + CONTROL_instance%EXCITE_SPECIES = "NONE" !$OMP PARALLEL private(nthreads, proc) proc = OMP_GET_THREAD_NUM() if(proc == 0) then @@ -1812,6 +1816,7 @@ subroutine CONTROL_load(unit) CONTROL_instance%DIIS_SWITCH_THRESHOLD_BKP = LowdinParameters_diisSwitchThreshold_bkp CONTROL_instance%ELECTRONIC_LEVEL_SHIFTING = LowdinParameters_electronicLevelShifting CONTROL_instance%NONELECTRONIC_LEVEL_SHIFTING = LowdinParameters_nonelectronicLevelShifting + CONTROL_instance%EXCHANGE_ORBITAL_THRESHOLD = LowdinParameters_exchangeOrbitalThreshold CONTROL_instance%WAVE_FUNCTION_SCALE = LowdinParameters_waveFunctionScale CONTROL_instance%SCF_NONELECTRONIC_MAX_ITERATIONS = LowdinParameters_scfNonelectronicMaxIterations CONTROL_instance%SCF_ELECTRONIC_MAX_ITERATIONS = LowdinParameters_scfElectronicMaxIterations @@ -2101,8 +2106,8 @@ subroutine CONTROL_load(unit) !! CONTROL_instance%MO_FRACTION_OCCUPATION = LowdinParameters_MOFractionOccupation CONTROL_instance%IONIZE_MO = LowdinParameters_ionizeMO - CONTROL_instance%IONIZE_SPECIE = LowdinParameters_ionizeSpecie - CONTROL_instance%EXCITE_SPECIE = LowdinParameters_exciteSpecie + CONTROL_instance%IONIZE_SPECIES = LowdinParameters_ionizeSpecies + CONTROL_instance%EXCITE_SPECIES = LowdinParameters_exciteSpecies !!***************************************************** !! Integrals transformation options @@ -2171,6 +2176,7 @@ subroutine CONTROL_save( unit, lastStep, firstStep ) LowdinParameters_diisSwitchThreshold_bkp = CONTROL_instance%DIIS_SWITCH_THRESHOLD_BKP LowdinParameters_electronicLevelShifting = CONTROL_instance%ELECTRONIC_LEVEL_SHIFTING LowdinParameters_nonelectronicLevelShifting = CONTROL_instance%NONELECTRONIC_LEVEL_SHIFTING + LowdinParameters_exchangeOrbitalThreshold = CONTROL_instance%EXCHANGE_ORBITAL_THRESHOLD LowdinParameters_waveFunctionScale = CONTROL_instance%WAVE_FUNCTION_SCALE LowdinParameters_scfNonelectronicMaxIterations = CONTROL_instance%SCF_NONELECTRONIC_MAX_ITERATIONS LowdinParameters_scfElectronicMaxIterations = CONTROL_instance%SCF_ELECTRONIC_MAX_ITERATIONS @@ -2449,8 +2455,8 @@ subroutine CONTROL_save( unit, lastStep, firstStep ) !! LowdinParameters_MOFractionOccupation = CONTROL_instance%MO_FRACTION_OCCUPATION LowdinParameters_ionizeMO = CONTROL_instance%IONIZE_MO - LowdinParameters_ionizeSpecie = CONTROL_instance%IONIZE_SPECIE - LowdinParameters_exciteSpecie = CONTROL_instance%EXCITE_SPECIE + LowdinParameters_ionizeSpecies = CONTROL_instance%IONIZE_SPECIES + LowdinParameters_exciteSpecies = CONTROL_instance%EXCITE_SPECIES !!***************************************************** !! Integrals transformation options @@ -2519,6 +2525,7 @@ subroutine CONTROL_copy(this, otherThis) otherThis%DIIS_SWITCH_THRESHOLD_BKP = this%DIIS_SWITCH_THRESHOLD_BKP otherThis%ELECTRONIC_LEVEL_SHIFTING = this%ELECTRONIC_LEVEL_SHIFTING otherThis%NONELECTRONIC_LEVEL_SHIFTING = this%NONELECTRONIC_LEVEL_SHIFTING + otherthis%EXCHANGE_ORBITAL_THRESHOLD = this%EXCHANGE_ORBITAL_THRESHOLD otherThis%WAVE_FUNCTION_SCALE= this%WAVE_FUNCTION_SCALE otherThis%SCF_NONELECTRONIC_MAX_ITERATIONS = this%SCF_NONELECTRONIC_MAX_ITERATIONS otherThis%SCF_ELECTRONIC_MAX_ITERATIONS = this%SCF_ELECTRONIC_MAX_ITERATIONS @@ -2764,8 +2771,8 @@ subroutine CONTROL_copy(this, otherThis) !! otherThis%MO_FRACTION_OCCUPATION = this%MO_FRACTION_OCCUPATION otherThis%IONIZE_MO = this%IONIZE_MO - otherThis%IONIZE_SPECIE = this%IONIZE_SPECIE - otherThis%EXCITE_SPECIE = this%EXCITE_SPECIE + otherThis%IONIZE_SPECIES = this%IONIZE_SPECIES + otherThis%EXCITE_SPECIES = this%EXCITE_SPECIES !!***************************************************** !! Integrals transformation options @@ -2795,7 +2802,7 @@ end subroutine CONTROL_copy subroutine CONTROL_show() implicit none - integer:: i !, nthreads, proc + integer:: i,j !, nthreads, proc print *,"" print *,"LOWDIN IS RUNNING WITH NEXT PARAMETERS: " @@ -2975,15 +2982,17 @@ subroutine CONTROL_show() end if - if((CONTROL_instance%IONIZE_SPECIE(1)) /= "NONE") then - - write (*,"(T10,A,I5)") "MOLECULAR ORBITAL TO BE IONIZED: ", CONTROL_instance%IONIZE_MO + if((CONTROL_instance%IONIZE_SPECIES(1)) /= "NONE") then ! print *, "size ionizepsecie", size(CONTROL_instance%IONIZE_SPECIE) ! print *, "ionizepsecie", CONTROL_instance%IONIZE_SPECIE - do i = 1, size(CONTROL_instance%IONIZE_SPECIE) - if ( trim(CONTROL_instance%IONIZE_SPECIE(i)) /= "NONE" ) then - write (*,"(T10,A,A)") "FOR SPECIE0: ", trim(CONTROL_instance%IONIZE_SPECIE(i)) - write (*,"(T10,A,ES15.5)") "IONIZED MOLECULAR ORBITAL OCCUPATION: ",CONTROL_instance%MO_FRACTION_OCCUPATION + write (*,"(T10,A)") "MOLECULAR ORBITALS TO BE IONIZED " + do i = 1, size(CONTROL_instance%IONIZE_SPECIES) + if ( trim(CONTROL_instance%IONIZE_SPECIES(i)) /= "NONE" ) then + write (*,"(T15,A,A)") "FOR SPECIES: ", trim(CONTROL_instance%IONIZE_SPECIES(i)) + do j = 1, size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(j) .gt. 0) & + write (*,"(T20,A,I4,A,ES15.5)") "ORBITAL", CONTROL_instance%IONIZE_MO(j) ," OCCUPATION:",CONTROL_instance%MO_FRACTION_OCCUPATION(j) + end do end if end do end if @@ -3000,9 +3009,9 @@ subroutine CONTROL_show() end if - if(CONTROL_instance%EXCITE_SPECIE /= "NONE") then + if(CONTROL_instance%EXCITE_SPECIES /= "NONE") then - write (*,"(T10,A,T10)") "CALCULATING", trim(CONTROL_instance%EXCITE_SPECIE) ,"IN THE FIRST EXCITED STATE" + write (*,"(T10,A,T10)") "CALCULATING", trim(CONTROL_instance%EXCITE_SPECIES) ,"IN THE FIRST EXCITED STATE" end if @@ -3157,6 +3166,12 @@ subroutine CONTROL_show() end if + if ( CONTROL_instance%EXCHANGE_ORBITALS_IN_SCF .eqv. .true. ) then + write(*,"(T10,A)") "OCCUPIED ORBITALS IN SCF WILL BE REORDERED TO MAXIMIZE THE OVERLAP WITH RESPECT TO THE INITIAL GUESS" + write(*,"(T10,A,F8.3)") "OVERLAP THRESHOLD FOR ORBITAL EXCHANGE:", CONTROL_instance%EXCHANGE_ORBITAL_THRESHOLD + end if + + if(CONTROL_instance%LOCALIZE_ORBITALS) then select case (trim(CONTROL_instance%ERKALE_LOCALIZATION_METHOD)) case("MU") diff --git a/src/core/Matrix.f90 b/src/core/Matrix.f90 index ef286ebe..b089a3b8 100644 --- a/src/core/Matrix.f90 +++ b/src/core/Matrix.f90 @@ -2659,16 +2659,108 @@ function Matrix_function( this, func, param ) result ( output ) end function Matrix_function + !> + !! @brief Calcula una funcion generica de la matriz this + function Matrix_function_svd( this, func, param ) result ( output ) + implicit none + type(Matrix), intent(inout) :: this + character(*), intent(in) :: func + type(Matrix) :: output + real(8), intent(in), optional :: param + + + type(Matrix) :: U, VT, singular + integer :: i + integer :: dim + + dim=Matrix_getNumberOfColumns(this) + + call Matrix_constructor(U, int(dim,8), int(dim,8), 0.0_8) + call Matrix_constructor(VT, int(dim,8), int(dim,8), 0.0_8) + call Matrix_constructor(singular, int(dim,8), int(dim,8), 0.0_8) + call Matrix_svd( this, U, VT, singular ) + + ! vectorsInverted = Matrix_inverse( range ) + + do i=1, size( singular%values, DIM=1 ) + + select case ( trim(func) ) + + case ( "pow" ) + singular%values(i, i) = singular%values(i, i)**param + + case ( "sqrt" ) + singular%values(i, i) = sqrt( singular%values(i, i) ) + + case ( "log" ) + singular%values(i, i) = log( singular%values(i, i) ) + + case ( "log10" ) + singular%values(i, i) = log10( singular%values(i, i) ) + + case ( "sin" ) + singular%values(i, i) = sin( singular%values(i, i) ) + + case ( "cos" ) + singular%values(i, i) = cos( singular%values(i, i) ) + + case ( "tan" ) + singular%values(i, i) = tan( singular%values(i, i) ) + + case ( "asin" ) + singular%values(i, i) = asin( singular%values(i, i) ) + + case ( "acos" ) + singular%values(i, i) = acos( singular%values(i, i) ) + + case ( "atan" ) + singular%values(i, i) = atan( singular%values(i, i) ) + + case ( "sinh" ) + singular%values(i, i) = sinh( singular%values(i, i) ) + + case ( "cosh" ) + singular%values(i, i) = cosh( singular%values(i, i) ) + + case ( "tanh" ) + singular%values(i, i) = tanh( singular%values(i, i) ) + + end select + + end do + + output = Matrix_product( singular, VT ) + output = Matrix_product( U, output ) + + end function Matrix_function_svd + !> !! @brief Calcula la potencia de la matriz this !! @warning La matriz tiene que ser diagonalizable - function Matrix_pow( this, eexponent ) result ( output ) + function Matrix_pow( this, eexponent, method ) result ( output ) implicit none type(Matrix), intent(inout) :: this real(8), intent(in) :: eexponent + character(*), optional :: method type(Matrix) :: output + character(10) :: selectedMethod + + if( present ( method ) ) then + selectedMethod=trim(method) + else + selectedMethod="eigen" + end if + + select case( trim(selectedMethod) ) + case("SVD") + output = Matrix_function_svd( this, "pow", eexponent ) - output = Matrix_function( this, "pow", eexponent ) + case("eigen") + output = Matrix_function( this, "pow", eexponent ) + + case default + call Matrix_exception(ERROR, "The selected method to compute the matrix power is not implemented", "Class object Matrix in the Matrix_pow() function") + end select end function Matrix_pow diff --git a/src/core/MolecularSystem.f90 b/src/core/MolecularSystem.f90 index 1aeff224..e141d5ad 100644 --- a/src/core/MolecularSystem.f90 +++ b/src/core/MolecularSystem.f90 @@ -198,9 +198,9 @@ subroutine MolecularSystem_build() end if !!Check for input errors in the number of particles - if( (abs(int(MolecularSystem_instance%species(i)%ocupationNumber)-MolecularSystem_instance%species(i)%ocupationNumber) & - .gt. CONTROL_instance%DOUBLE_ZERO_THRESHOLD) .and. (CONTROL_instance%MO_FRACTION_OCCUPATION .eq. 1.0_8) ) then - print *, "species ", trim(MolecularSystem_getNameOfSpecie(i)) , "has fractional ocupation number ", MolecularSystem_instance%species(i)%ocupationNumber, "please check your input addParticles and multiplicity" + if( (abs(int(MolecularSystem_instance%species(i)%ocupationNumber)-MolecularSystem_instance%species(i)%ocupationNumber) .gt. CONTROL_instance%DOUBLE_ZERO_THRESHOLD)) then + print *, "species ", trim(MolecularSystem_getNameOfSpecie(i)) , "has fractional ocupation number ", & + MolecularSystem_instance%species(i)%ocupationNumber, "please check your input addParticles and multiplicity" call MolecularSystem_exception(ERROR, "Fractional ocupation number, imposible combination of charge and multiplicity","MolecularSystem module at build function.") end if end do @@ -333,15 +333,23 @@ end subroutine MolecularSystem_rotateOnPrincipalAxes !> !! @brief shows general information of the molecular system. !! @author S. A. Gonzalez - subroutine MolecularSystem_showInformation() - implicit none + subroutine MolecularSystem_showInformation(this) + type(MolecularSystem), optional, target :: this + + type(MolecularSystem), pointer :: system + + if( present(this) ) then + system=>this + else + system=>MolecularSystem_instance + end if print *,"" - print *," MOLECULAR SYSTEM: ",trim(MolecularSystem_instance%name) + print *," MOLECULAR SYSTEM: ",trim(system%name) print *,"-----------------" print *,"" - write (6,"(T5,A16,A)") "DESCRIPTION : ", trim( MolecularSystem_instance%description ) - write (6,"(T5,A16,I3)") "CHARGE : ",MolecularSystem_instance%charge + 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" print *,"" @@ -353,18 +361,26 @@ end subroutine MolecularSystem_showInformation !! @author S. A. Gonzalez !! @par changes : !! - rewritten, E. F. Posada. 2013 - subroutine MolecularSystem_showParticlesInformation() + subroutine MolecularSystem_showParticlesInformation(this) implicit none - integer :: i - integer :: j + type(MolecularSystem), optional, target :: this + type(MolecularSystem), pointer :: system + integer :: i, j + + if( present(this) ) then + system=>this + else + system=>MolecularSystem_instance + end if + print *,"" print *," INFORMATION OF PARTICLES :" print *,"=========================" - write (6,"(T10,A29,I8.0)") "Total number of particles =", MolecularSystem_instance%numberOfQuantumParticles + MolecularSystem_instance%numberOfPointCharges - write (6,"(T10,A29,I8.0)") "Number of quantum particles =", MolecularSystem_instance%numberOfQuantumParticles - write (6,"(T10,A29,I8.0)") "Number of puntual charges =", MolecularSystem_instance%numberOfPointCharges - write (6,"(T10,A29,I8.0)") "Number of quantum species =", MolecularSystem_instance%numberOfQuantumSpecies + write (6,"(T10,A29,I8.0)") "Total number of particles =", system%numberOfQuantumParticles + system%numberOfPointCharges + write (6,"(T10,A29,I8.0)") "Number of quantum particles =", system%numberOfQuantumParticles + write (6,"(T10,A29,I8.0)") "Number of puntual charges =", system%numberOfPointCharges + write (6,"(T10,A29,I8.0)") "Number of quantum species =", system%numberOfQuantumSpecies !!*********************************************************************** !! Imprime iformacion sobre masa, carga y numero de particulas encontradas @@ -375,29 +391,29 @@ subroutine MolecularSystem_showParticlesInformation() write (6,"(T10,A2,A4,A8,A12,A4,A5,A6,A5,A4,A5,A12)") "ID", " ","Symbol", " ","mass", " ","charge", " ","spin","","multiplicity" write (6,"(T5,A70)") "---------------------------------------------------------------------" - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + do i = 1, system%numberOfQuantumSpecies write (6,'(T8,I3.0,A5,A10,A5,F10.4,A5,F5.2,A5,F5.2,A5,F5.2)') & i, " ", & - trim(MolecularSystem_instance%species(i)%symbol)," ",& - MolecularSystem_instance%species(i)%mass," ",& - MolecularSystem_instance%species(i)%charge, " ",& - MolecularSystem_instance%species(i)%spin, "",& - MolecularSystem_instance%species(i)%multiplicity + trim(system%species(i)%symbol)," ",& + system%species(i)%mass," ",& + system%species(i)%charge, " ",& + system%species(i)%spin, "",& + system%species(i)%multiplicity end do print *,"" print *," CONSTANTS OF COUPLING " write (6,"(T7,A60)") "------------------------------------------------------------" - write (6,"(T10,A8,A10,A5,A5,A4,A5,A6,A5,A10)") "Symbol", " ","kappa", " ","eta", " ","lambda","","occupation" + write (6,"(T10,A11,A11,A11,A11,A11)") "Symbol", "kappa","eta","lambda","occupation" write (6,"(T7,A60)") "------------------------------------------------------------" - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies - write (6,'(T10,A10,A5,F7.1,A5,F5.2,A5,F5.2,A5,F5.2,A5,F5.2)') & - trim(MolecularSystem_instance%species(i)%symbol)," ",& - MolecularSystem_instance%species(i)%kappa, " ", & - MolecularSystem_instance%species(i)%eta, " ",& - MolecularSystem_instance%species(i)%lambda," ",& - MolecularSystem_instance%species(i)%ocupationNumber + do i = 1, system%numberOfQuantumSpecies + write (6,'(T10,A11,F11.2,F11.2,F11.2,F11.2)') & + trim(system%species(i)%symbol),& + system%species(i)%kappa,& + system%species(i)%eta,& + system%species(i)%lambda,& + system%species(i)%ocupationNumber end do print *,"" @@ -408,26 +424,26 @@ subroutine MolecularSystem_showParticlesInformation() write (6,"(T7,A60)") "------------------------------------------------------------" !! Only shows the basis-set name of the first particle by specie. - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + do i = 1, system%numberOfQuantumSpecies - if( MolecularSystem_instance%species(i)%isElectron .and. CONTROL_instance%IS_OPEN_SHELL ) then + if( system%species(i)%isElectron .and. CONTROL_instance%IS_OPEN_SHELL ) then write (6,'(T10,A10,A5,I8,A5,I12,A5,A10)') & - trim(MolecularSystem_instance%species(i)%symbol)," ",& + trim(system%species(i)%symbol)," ",& !MolecularSystem_getTotalNumberOfContractions(i)," ",& - MolecularSystem_instance%species(i)%basisSetSize," ",& - int(MolecularSystem_instance%species(i)%internalSize / 2), " ",& - trim(MolecularSystem_instance%species(i)%particles(1)%basis%name) + system%species(i)%basisSetSize," ",& + int(system%species(i)%internalSize / 2), " ",& + trim(system%species(i)%particles(1)%basis%name) ! write(*,*)MolecularSystem_getTotalNumberOfContractions(i) else write (6,'(T10,A10,A5,I8,A5,I12,A5,A10)') & - trim(MolecularSystem_instance%species(i)%symbol)," ",& + trim(system%species(i)%symbol)," ",& !MolecularSystem_getTotalNumberOfContractions(i)," ",& - MolecularSystem_instance%species(i)%basisSetSize," ",& - MolecularSystem_instance%species(i)%internalSize, " ",& - trim(MolecularSystem_instance%species(i)%particles(1)%basis%name) + system%species(i)%basisSetSize," ",& + system%species(i)%internalSize, " ",& + trim(system%species(i)%particles(1)%basis%name) ! write(*,*)MolecularSystem_getTotalNumberOfContractions(i) end if @@ -440,20 +456,20 @@ subroutine MolecularSystem_showParticlesInformation() write (6,"(T10,A11,A9,A20,A20)") " PRIMITIVE ", " SHELL "," EXPONENT "," COEFFICIENT " write (6,"(T10,A60)") "------------------------------------------------------------" - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + do i = 1, system%numberOfQuantumSpecies !! Avoid print twice basis in open-shell case - if(trim(MolecularSystem_instance%species(i)%name) == "E-BETA" ) cycle + if(trim(system%species(i)%name) == "E-BETA" ) cycle write (6,*) "" - write( 6, "(T5,A32,A5)") "BEGIN DESCRIPTION OF BASIS FOR: ", trim(MolecularSystem_instance%species(i)%symbol) + write( 6, "(T5,A32,A5)") "BEGIN DESCRIPTION OF BASIS FOR: ", trim(system%species(i)%symbol) write (6,"(T5,A30)") "================================" write (6,*) "" - do j = 1, size( MolecularSystem_instance%species(i)%particles ) + do j = 1, size( system%species(i)%particles ) - call BasisSet_showInCompactForm( MolecularSystem_instance%species(i)%particles(j)%basis,& - trim(MolecularSystem_instance%species(i)%particles(j)%nickname )) + call BasisSet_showInCompactForm( system%species(i)%particles(j)%basis,& + trim(system%species(i)%particles(j)%nickname )) end do @@ -474,9 +490,9 @@ subroutine MolecularSystem_showParticlesInformation() write (6,"(T10,A10,A5,A17,A5,A10)") "ID", " ", "Occupied Orbitals", " " ,"Basis size" write (6,"(T5,A70)") "---------------------------------------------------------------------" - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + do i = 1, system%numberOfQuantumSpecies write (6,'(T10,A10,A5,I8,A5,I12)') & - trim(MolecularSystem_instance%species(i)%symbol)," ",& + trim(system%species(i)%symbol)," ",& MolecularSystem_getOcupationNumber( i )," ",& MolecularSystem_getTotalNumberOfContractions(i) end do @@ -599,18 +615,24 @@ end subroutine MolecularSystem_showDistanceMatrix !> !! @brief Saves all system information. !! @author E. F. Posada, 2013 - subroutine MolecularSystem_saveToFile() + subroutine MolecularSystem_saveToFile(targetFilePrefix) implicit none - integer i, j, k + character(*), optional :: targetFilePrefix + character(50) :: fileName, prefix character(100) :: title + integer i, j, k + + prefix="lowdin" + if ( present( targetFilePrefix ) ) prefix=trim(targetFilePrefix) !!**************************************************************************** !! CONTROL parameters on file. !! !! open file - open(unit=40, file="lowdin.dat", status="replace", form="formatted") + filename=trim(prefix)//".dat" + open(unit=40, file=filename, status="replace", form="formatted") !!save all options call CONTROL_save(40) @@ -621,7 +643,8 @@ subroutine MolecularSystem_saveToFile() !! !!Open file - open(unit=40, file="lowdin.sys", status="replace", form="formatted") + filename=trim(prefix)//".sys" + open(unit=40, file=filename, status="replace", form="formatted") !! Saving general information. write(40,*) MolecularSystem_instance%name @@ -674,7 +697,8 @@ subroutine MolecularSystem_saveToFile() !! !!open file - open(unit=40, file="lowdin.bas", status="replace", form="formatted") + filename=trim(prefix)//".bas" + open(unit=40, file=filename, status="replace", form="formatted") write(40,*) MolecularSystem_instance%numberOfQuantumSpecies do i = 1, MolecularSystem_instance%numberOfQuantumSpecies @@ -743,12 +767,12 @@ end subroutine MolecularSystem_saveToFile !> !! @brief loads all system information from file !! @author E. F. Posada, 2013 - subroutine MolecularSystem_loadFromFile( form, targetFile ) + subroutine MolecularSystem_loadFromFile( form, targetPrefix ) implicit none character(*) :: form - character(*), optional :: targetFile - character(50) :: fileName + character(*), optional :: targetPrefix + character(50) :: filePrefix,fileName integer :: auxValue integer :: counter @@ -758,91 +782,86 @@ subroutine MolecularSystem_loadFromFile( form, targetFile ) character(50) :: species character(50) :: otherSpecies - fileName="lowdin" - if ( present( targetFile ) ) fileName=trim(targetFile) + filePrefix="lowdin" + if ( present( targetPrefix ) ) filePrefix=trim(targetPrefix) + + existFile = .false. select case (trim(form)) case("LOWDIN.BAS") - + !!**************************************************************************** !! Loading info from the lowdin.bas format !! - + !!open file - existFile = .false. - inquire(file="lowdin.bas", exist=existFile) - - if(existFile) then + fileName=trim(filePrefix)//".bas" + inquire(file=fileName, exist=existFile) + if( .not. existFile) call MolecularSystem_exception(ERROR, "The file: "//trim(fileName)//" was not found!","MolecularSystem module at LoadFromFile function.") - !! Destroy the molecular system if any - ! call MolecularSystem_destroy() - if(allocated(MolecularSystem_instance%pointCharges)) deallocate(MolecularSystem_instance%pointCharges) - if(allocated(MolecularSystem_instance%allParticles)) deallocate(MolecularSystem_instance%allParticles) + !! Destroy the molecular system if any + ! call MolecularSystem_destroy() + if(allocated(MolecularSystem_instance%pointCharges)) deallocate(MolecularSystem_instance%pointCharges) + if(allocated(MolecularSystem_instance%allParticles)) deallocate(MolecularSystem_instance%allParticles) - open(unit=40, file="lowdin.bas", status="old", form="formatted") - - read(40,*) MolecularSystem_instance%numberOfQuantumSpecies - if(.not. allocated(MolecularSystem_instance%species)) allocate(MolecularSystem_instance%species(MolecularSystem_instance%numberOfQuantumSpecies)) + open(unit=40, file=fileName, status="old", form="formatted") - MolecularSystem_instance%numberOfQuantumParticles = 0 + read(40,*) MolecularSystem_instance%numberOfQuantumSpecies + if(.not. allocated(MolecularSystem_instance%species)) allocate(MolecularSystem_instance%species(MolecularSystem_instance%numberOfQuantumSpecies)) - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies - - if(allocated(MolecularSystem_instance%species(i)%particles)) deallocate(MolecularSystem_instance%species(i)%particles) + MolecularSystem_instance%numberOfQuantumParticles = 0 - read(40,*) MolecularSystem_instance%species(i)%name - read(40,*) auxValue - - allocate(MolecularSystem_instance%species(i)%particles(auxValue)) - - do j = 1, size(MolecularSystem_instance%species(i)%particles) - - read(40,*) MolecularSystem_instance%species(i)%particles(j)%nickname - - MolecularSystem_instance%numberOfQuantumParticles = MolecularSystem_instance%numberOfQuantumParticles + 1 - call BasisSet_load(MolecularSystem_instance%species(i)%particles(j)%basis, "LOWDIN.BAS", unit = 40) - - end do - - MolecularSystem_instance%species(i)%basisSetSize = MolecularSystem_getNumberOfContractions(i) - - end do + do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + + if(allocated(MolecularSystem_instance%species(i)%particles)) deallocate(MolecularSystem_instance%species(i)%particles) + + read(40,*) MolecularSystem_instance%species(i)%name + read(40,*) auxValue + + allocate(MolecularSystem_instance%species(i)%particles(auxValue)) + + do j = 1, size(MolecularSystem_instance%species(i)%particles) + + read(40,*) MolecularSystem_instance%species(i)%particles(j)%nickname + + MolecularSystem_instance%numberOfQuantumParticles = MolecularSystem_instance%numberOfQuantumParticles + 1 + call BasisSet_load(MolecularSystem_instance%species(i)%particles(j)%basis, filename, unit = 40) - read(40,*) MolecularSystem_instance%numberOfPointCharges - allocate(MolecularSystem_instance%pointCharges(MolecularSystem_instance%numberOfPointCharges)) - - do i = 1, MolecularSystem_instance%numberOfPointCharges - read(40,*) MolecularSystem_instance%pointCharges(i)%charge - read(40,*) MolecularSystem_instance%pointCharges(i)%origin - end do - - close(40) - - !! Set the particles manager (all pointers) - MolecularSystem_instance%numberOfParticles = MolecularSystem_instance%numberOfQuantumParticles + MolecularSystem_instance%numberOfPointCharges - allocate(molecularSystem_instance%allParticles(MolecularSystem_instance%numberOfParticles )) - - counter = 1 - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies - do j = 1, size(MolecularSystem_instance%species(i)%particles) - molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%species(i)%particles(j) - counter = counter + 1 - end do end do - do i = 1, MolecularSystem_instance%numberOfPointCharges - molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%pointCharges(i) + MolecularSystem_instance%species(i)%basisSetSize = MolecularSystem_getNumberOfContractions(i) + + end do + + read(40,*) MolecularSystem_instance%numberOfPointCharges + allocate(MolecularSystem_instance%pointCharges(MolecularSystem_instance%numberOfPointCharges)) + + do i = 1, MolecularSystem_instance%numberOfPointCharges + read(40,*) MolecularSystem_instance%pointCharges(i)%charge + read(40,*) MolecularSystem_instance%pointCharges(i)%origin + end do + + close(40) + + !! Set the particles manager (all pointers) + MolecularSystem_instance%numberOfParticles = MolecularSystem_instance%numberOfQuantumParticles + MolecularSystem_instance%numberOfPointCharges + allocate(molecularSystem_instance%allParticles(MolecularSystem_instance%numberOfParticles )) + + counter = 1 + do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + do j = 1, size(MolecularSystem_instance%species(i)%particles) + molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%species(i)%particles(j) counter = counter + 1 end do - - particleManager_instance => molecularSystem_instance%allParticles + end do - else - - call MolecularSystem_exception(ERROR, "The file: "//trim(fileName)//".bas was not found!","MolecularSystem module at LoadFromFile function.") - - end if + do i = 1, MolecularSystem_instance%numberOfPointCharges + molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%pointCharges(i) + counter = counter + 1 + end do + + particleManager_instance => molecularSystem_instance%allParticles case("LOWDIN.DAT") @@ -851,24 +870,17 @@ subroutine MolecularSystem_loadFromFile( form, targetFile ) !! !! open file - existFile = .false. - inquire(file="lowdin.dat", exist=existFile) + fileName=trim(filePrefix)//".dat" + inquire(file=fileName, exist=existFile) + if( .not. existFile) call MolecularSystem_exception(ERROR, "The file: "//trim(fileName)//" was not found!","MolecularSystem module at LoadFromFile function.") - if(existFile) then - - open(unit=40, file="lowdin.dat", status="old", form="formatted") - - call CONTROL_start() - call CONTROL_load(unit = 40) - - close(40) - - else - - call MolecularSystem_exception(ERROR, "The file: "//trim(fileName)//".dat was not found!","MolecularSystem module at LoadFromFile function.") - - end if + open(unit=40, file=fileName, status="old", form="formatted") + call CONTROL_start() + call CONTROL_load(unit = 40) + + close(40) + case("LOWDIN.SYS") !!**************************************************************************** @@ -879,107 +891,101 @@ subroutine MolecularSystem_loadFromFile( form, targetFile ) call MolecularSystem_destroy() !! open file - existFile = .false. - inquire(file="lowdin.sys", exist=existFile) - - if(existFile) then - !!Open file - open(unit=40, file="lowdin.sys", status="old", form="formatted") - - !! read general information. - read(40,*) MolecularSystem_instance%name - read(40,'(A100)') MolecularSystem_instance%description - read(40,*) MolecularSystem_instance%charge - - !! load quantum species. - read(40,*) MolecularSystem_instance%numberOfQuantumSpecies - allocate(MolecularSystem_instance%species(MolecularSystem_instance%numberOfQuantumSpecies)) + fileName=trim(filePrefix)//".sys" + inquire(file=fileName, exist=existFile) + if( .not. existFile) call MolecularSystem_exception(ERROR, "The file: "//trim(fileName)//" was not found!","MolecularSystem module at LoadFromFile function.") - !! load particles for each species. - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies - call Species_loadFromFile(MolecularSystem_instance%species(i), unit=40) - end do + !!Open file + open(unit=40, file=fileName, status="old", form="formatted") - !! load Point charges - read(40,*) MolecularSystem_instance%numberOfPointCharges - allocate(MolecularSystem_instance%pointCharges(MolecularSystem_instance%numberOfPointCharges)) - - !! load info of each point charge - do i = 1, MolecularSystem_instance%numberOfPointCharges - call Particle_loadFromFile(MolecularSystem_instance%pointCharges(i), unit=40) - end do - - !! read the total of particles on the system - read(40,*) MolecularSystem_instance%numberOfParticles - read(40,*) MolecularSystem_instance%numberOfQuantumParticles - + !! read general information. + read(40,*) MolecularSystem_instance%name + read(40,'(A100)') MolecularSystem_instance%description + read(40,*) MolecularSystem_instance%charge - !! Set the particles manager (all pointers) - allocate(molecularSystem_instance%allParticles(MolecularSystem_instance%numberOfParticles )) - - counter = 1 - do i = 1, MolecularSystem_instance%numberOfQuantumSpecies - do j = 1, size(MolecularSystem_instance%species(i)%particles) - - molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%species(i)%particles(j) - counter = counter + 1 - - end do - end do - - do i = 1, MolecularSystem_instance%numberOfPointCharges - - molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%pointCharges(i) + !! load quantum species. + read(40,*) MolecularSystem_instance%numberOfQuantumSpecies + allocate(MolecularSystem_instance%species(MolecularSystem_instance%numberOfQuantumSpecies)) + + !! load particles for each species. + do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + call Species_loadFromFile(MolecularSystem_instance%species(i), unit=40) + end do + + !! load Point charges + read(40,*) MolecularSystem_instance%numberOfPointCharges + allocate(MolecularSystem_instance%pointCharges(MolecularSystem_instance%numberOfPointCharges)) + + !! load info of each point charge + do i = 1, MolecularSystem_instance%numberOfPointCharges + call Particle_loadFromFile(MolecularSystem_instance%pointCharges(i), unit=40) + end do + + !! read the total of particles on the system + read(40,*) MolecularSystem_instance%numberOfParticles + read(40,*) MolecularSystem_instance%numberOfQuantumParticles + + + !! Set the particles manager (all pointers) + allocate(molecularSystem_instance%allParticles(MolecularSystem_instance%numberOfParticles )) + + counter = 1 + do i = 1, MolecularSystem_instance%numberOfQuantumSpecies + do j = 1, size(MolecularSystem_instance%species(i)%particles) + + molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%species(i)%particles(j) counter = counter + 1 - + end do - - particleManager_instance => molecularSystem_instance%allParticles + end do - - !! Loading External/Inter-particle potentials information - if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) then + do i = 1, MolecularSystem_instance%numberOfPointCharges + + molecularSystem_instance%allParticles(counter)%particlePtr => MolecularSystem_instance%pointCharges(i) + counter = counter + 1 - read(40,*) auxValue - call ExternalPotential_constructor(auxValue) + end do - !! FELIX TODO: create function to get potential ID - - do j = 1, ExternalPotential_instance%ssize - read(40,*) i - read(40,*) name - read(40,*) species + particleManager_instance => molecularSystem_instance%allParticles - call ExternalPotential_load(i, name, species) - end do + !! Loading External/Inter-particle potentials information + if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) then - end if + read(40,*) auxValue + call ExternalPotential_constructor(auxValue) - if(CONTROL_instance%IS_THERE_INTERPARTICLE_POTENTIAL) then + !! FELIX TODO: create function to get potential ID - read(40,*) auxValue - call InterPotential_constructor(auxValue) + do j = 1, ExternalPotential_instance%ssize + read(40,*) i + read(40,*) name + read(40,*) species - do j = 1, InterPotential_instance%ssize - read(40,*) i - read(40,*) name - read(40,*) species - read(40,*) otherSpecies - - call InterPotential_load(i, name, species, otherSpecies) + call ExternalPotential_load(i, name, species) - end do + end do - end if + end if - close(40) + if(CONTROL_instance%IS_THERE_INTERPARTICLE_POTENTIAL) then + + read(40,*) auxValue + call InterPotential_constructor(auxValue) + + do j = 1, InterPotential_instance%ssize + read(40,*) i + read(40,*) name + read(40,*) species + read(40,*) otherSpecies + + call InterPotential_load(i, name, species, otherSpecies) + + end do - else - - call MolecularSystem_exception(ERROR, "The file: "//trim(fileName)//".sys was not found!","MolecularSystem module at LoadFromFile function.") - end if + + close(40) end select diff --git a/src/core/Vector.f90 b/src/core/Vector.f90 index 35796f79..b199c8a9 100644 --- a/src/core/Vector.f90 +++ b/src/core/Vector.f90 @@ -110,7 +110,9 @@ module Vector_ Vector_constructorInteger, & Vector_constructorInteger8, & Vector_destructorInteger, & - Vector_swapIntegerElements + Vector_swapIntegerElements, & + Vector_writeToFileInteger, & + Vector_getFromFileInteger contains @@ -313,7 +315,7 @@ end subroutine Vector_constructorInteger8 !> !! @brief Constructor de copia - !! Reserva la memoria necesaria para otherMatrix y le asigna los valores de this + !! Reserva la memoria necesaria para otherVector y le asigna los valores de this subroutine Vector_copyConstructorInteger( this, otherVector ) implicit none type(IVector), intent(inout) :: this @@ -330,7 +332,7 @@ end subroutine Vector_copyConstructorInteger !> !! @brief Constructor de copia - !! Reserva la memoria necesaria para otherMatrix y le asigna los valores de this + !! Reserva la memoria necesaria para otherVector y le asigna los valores de this subroutine Vector_copyConstructor( this, otherVector ) implicit none type(Vector), intent(inout) :: this @@ -345,7 +347,7 @@ end subroutine Vector_copyConstructor !> !! @brief Constructor de copia - !! Reserva la memoria necesaria para otherMatrix y le asigna los valores de this + !! Reserva la memoria necesaria para otherVector y le asigna los valores de this subroutine Vector_copyConstructor8( this, otherVector ) implicit none type(Vector8), intent(inout) :: this @@ -403,7 +405,6 @@ subroutine Vector_show( this, flags, keys ) integer :: columns integer :: ssize integer :: i - integer :: j integer :: k integer :: l integer :: u @@ -500,7 +501,6 @@ subroutine Vector_showInteger( this, flags, keys ) integer :: columns integer :: ssize integer :: i - integer :: j integer :: k integer :: l integer :: u @@ -599,7 +599,6 @@ subroutine Vector_writeToFile(vvector, unit, file, binary, value, arguments) character(*), optional :: arguments(:) integer :: elementsNum - integer :: status integer :: n character(20) :: auxSize @@ -640,7 +639,7 @@ subroutine Vector_writeToFile(vvector, unit, file, binary, value, arguments) else call Vector_exception( ERROR, "Unit file no connected!",& - "Class object Matrix in the writeToFile() function" ) + "Class object Vector in the writeToFile() function" ) end if @@ -694,7 +693,7 @@ subroutine Vector_writeToFile(vvector, unit, file, binary, value, arguments) else call Vector_exception( ERROR, "Unit file no connected!",& - "Class object Matrix in the writeToFile() function" ) + "Class object Vector in the writeToFile() function" ) end if @@ -717,7 +716,6 @@ subroutine Vector_getFromFile(elementsNum, unit, file, binary, value, arguments, real(8), optional :: value type(Vector), optional, intent(out) :: output - type(Exception) :: ex character(5000) :: line character(20) :: auxSize integer :: status @@ -746,11 +744,11 @@ subroutine Vector_getFromFile(elementsNum, unit, file, binary, value, arguments, rewind(unit) found = .false. - line = "" if(present(arguments)) then do + line = "" read(unit, iostat = status) line (1:len_trim(arguments(1))) if(status == -1) then @@ -814,8 +812,8 @@ subroutine Vector_getFromFile(elementsNum, unit, file, binary, value, arguments, else - call Vector_exception( ERROR, "The dimensions of the matrix "//trim(file)//" are wrong ",& - "Class object Matrix in the getFromFile() function" ) + call Vector_exception( ERROR, "The dimensions of the vector "//trim(file)//" are wrong ",& + "Class object Vector in the getFromFile() function" ) end if @@ -824,7 +822,7 @@ subroutine Vector_getFromFile(elementsNum, unit, file, binary, value, arguments, else call Vector_exception( ERROR, "Unit file no connected!",& - "Class object Matrix in the getFromFile() function" ) + "Class object Vector in the getFromFile() function" ) end if @@ -962,8 +960,8 @@ subroutine Vector_getFromFile(elementsNum, unit, file, binary, value, arguments, else - call Vector_exception( ERROR, "The dimensions of the matrix "//trim(file)//" are wrong ",& - "Class object Matrix in the getFromFile() function" ) + call Vector_exception( ERROR, "The dimensions of the vector "//trim(file)//" are wrong ",& + "Class object Vector in the getFromFile() function" ) end if @@ -972,17 +970,407 @@ subroutine Vector_getFromFile(elementsNum, unit, file, binary, value, arguments, else call Vector_exception( ERROR, "Unit file no connected!",& - "Class object Matrix in the getFromFile() function" ) + "Class object Vector in the getFromFile() function" ) end if end if end if end subroutine Vector_getFromFile + + !> + !! @brief Escribe un vector en el lugar especificado (en binario o con formato) + subroutine Vector_writeToFileInteger(intvector, unit, file, binary, value, arguments) + implicit none + type(IVector), optional, intent(in) :: intvector + integer, optional :: unit + character(*), optional :: file + logical, optional :: binary + integer, optional :: value + character(*), optional :: arguments(:) + + integer :: elementsNum + integer :: n + + character(20) :: auxSize + + logical :: bbinary + logical :: existFile + + bbinary = .false. + if(present(binary)) bbinary = binary + + if ( bbinary ) then + if ( present( unit ) ) then + !! It is assumed that the unit y conected to any file (anyways will check) + inquire(unit=unit, exist=existFile) + + if(existFile) then + + if(present(arguments)) then + do n = 1, size(arguments) + write(unit) arguments(n) + end do + end if + + if(present(value)) then + write(unit) 1_8 + write(unit) value + + else + + write(unit) int(size(intvector%values), 8) + write(unit) intvector%values + + end if + + else + + call Vector_exception( ERROR, "Unit file no connected!",& + "Class object Vector in the writeToFile() function" ) + + end if + + + else if ( present(file) ) then + if(bbinary) then + open ( 4,FILE=trim(file),STATUS='REPLACE',ACTION='WRITE', FORM ='UNFORMATTED') + write(4) int(size(intvector%values), 8) + write(4) intvector%values + close(4) + + else + + open ( 4,FILE=trim(file),STATUS='REPLACE',ACTION='WRITE') + elementsNum = size( intvector%values ) + write(auxSize,*) elementsNum + write (4,"("//trim(auxSize)//"I15)") (intvector%values(n), n=1 , elementsNum) + close(4) + end if + + end if + + else !! not binary + + if ( present( unit ) ) then + !! It is assumed that the unit y conected to any file (anyways will check) + inquire(unit=unit, exist=existFile) + + if(existFile) then + + if(present(arguments)) then + + do n = 1, size(arguments) + + write(unit,*) arguments(n) + + end do + end if + + if(present(value)) then + write(unit,*) 1_8 + write(unit,*) value + + else + + write(unit,*) int(size(intvector%values), 8) + write(unit,*) intvector%values + + end if + + else + + call Vector_exception( ERROR, "Unit file no connected!",& + "Class object Vector in the writeToFile() function" ) + + end if + + end if !! unit not present + + + end if + + end subroutine Vector_writeToFileInteger + + !> + !! @brief Obtiene un vector del lugar especificado + subroutine Vector_getFromFileInteger(elementsNum, unit, file, binary, value, arguments, output ) + implicit none + integer, optional, intent(in) :: elementsNum + integer, optional :: unit + character(*), optional :: file + logical, optional :: binary + character(*), optional :: arguments(:) + integer, optional :: value + type(IVector), optional, intent(out) :: output + + character(5000) :: line + character(20) :: auxSize + integer :: status + integer :: n + integer(8) :: totalSize + logical :: bbinary + logical :: existFile + logical :: found + + + if (present(elementsNum)) write(auxSize,*) elementsNum + + bbinary = .false. + existFile = .false. + + if(present(binary)) bbinary = binary + + if ( bbinary ) then + if ( present( unit ) ) then + + !! check file + inquire(unit=unit, exist=existFile) + + if(existFile) then + + rewind(unit) + + found = .false. + line = "" + + if(present(arguments)) then + + do + read(unit, iostat = status) line (1:len_trim(arguments(1))) + + if(status == -1) then + + print *, "while searching for", arguments + call vector_exception( ERROR, "End of file!",& + "Class object Vector in the getfromFile() function") + end if + + if(trim(line) == trim(arguments(1))) then + + found = .true. + + end if + + if(found) then + + backspace(unit) + + do n = 1, size(arguments) + + found = .false. + read(unit, iostat = status) line + + if(trim(line) == trim(arguments(n))) then + + found = .true. + + end if + + end do + + end if + + if(found) exit + + end do + + + end if + + !! check size + read(unit) totalSize + + if(present(value)) then + + read(unit) value + + else + + if(totalSize == int(elementsNum,8)) then + + if(.not. allocated(output%values)) then + + call Vector_constructorInteger( output, elementsNum ) + + end if + + read(unit) output%values + + ! call Vector_show(output) + + else + + print *, "while searching for", arguments, "totalSize was", totalSize, "and expected", int(elementsNum,8) + call Vector_exception( ERROR, "The dimensions of the vector "//trim(file)//" are wrong ",& + "Class object Vector in the getFromFile() function" ) + + end if + + end if + + else + + call Vector_exception( ERROR, "Unit file no connected!",& + "Class object Vector in the getFromFile() function" ) + + end if + + else if ( present(file) ) then + + inquire( FILE = trim(file), EXIST = existFile ) + + if ( existFile ) then + + if(bbinary) then + open( 4, FILE=trim(file), ACTION='read', FORM='unformatted', STATUS='old' ) + + read(4) totalSize + if(totalSize == int(elementsNum,8)) then + + call Vector_constructorInteger( output, elementsNum ) + + read(4) output%values + close(4) + + !! print*, output%values + + else + close(4) + + call Vector_exception(ERROR, "The dimensions of the vector in "//trim(file)//" are wrong ", & + "Class object Vector_ in the getFromFile() function") + end if + else + + call Vector_constructorInteger( output, elementsNum ) + open ( 4,FILE=trim(file),STATUS='unknown',ACCESS='sequential' ) + !! verifica el tamahno de las matrice + read (4,"(A)") line + if ( len_trim(line) /=elementsNum*15) then + close(4) + + !! call tracebackqq(string="Probando trace ",USER_EXIT_CODE=-1,EPTR=%LOC(ExceptionInfo)) + + call Vector_exception(ERROR, "The dimensions of the vector in "//trim(file)//" are wrong ", & + "Class object Vector_ in the getFromFile() function" ) + + else + rewind(4) + end if + read ( 4,"("//trim(auxSize)//"I15)") (output%values(n), n=1 ,elementsNum) + close(4) + + end if + else + + call Vector_exception(ERROR, "The file "//trim(file)//" don't exist " , & + "Class object Vector_ in the getFromFile() function" ) + + end if + end if + + else !! not binary + + if ( present( unit ) ) then + + !! check file + inquire(unit=unit, exist=existFile) + + if(existFile) then + + rewind(unit) + + found = .false. + line = "" + + if(present(arguments)) then + + do + read(unit, *, iostat = status) line (1:len_trim(arguments(1))) + + if(status == -1) then + + call vector_exception( ERROR, "End of file!",& + "Class object Vector in the getfromFile() function" ) + end if + + if(trim(line) == trim(arguments(1))) then + + found = .true. + + end if + + if(found) then + + backspace(unit) + + do n = 1, size(arguments) + + found = .false. + read(unit, *, iostat = status) line + + if(trim(line) == trim(arguments(n))) then + + found = .true. + + end if + + end do + + end if + + if(found) exit + + end do + + + end if + + !! check size + read(unit,*) totalSize + + if(present(value)) then + + read(unit,*) value + + else + + if(totalSize == int(elementsNum,8)) then + + if(.not. allocated(output%values)) then + + call Vector_constructorInteger( output, elementsNum ) + + end if + + read(unit,*) output%values + + ! call Vector_show(output) + + else + + call Vector_exception( ERROR, "The dimensions of the vector "//trim(file)//" are wrong ",& + "Class object Vector in the getFromFile() function" ) + + end if + + end if + + else + + call Vector_exception( ERROR, "Unit file no connected!",& + "Class object Vector in the getFromFile() function" ) + + end if + end if + end if + + end subroutine Vector_getFromFileInteger !> - !! @brief Devuelve un apuntador a la matrix solicitada - !! @param this matrix de m x n + !! @brief Devuelve un apuntador a la vector solicitada + !! @param this vector de m x n !! @return Apuntador a la matriz solicitada. !! @todo No ha sido probada function Vector_getPtr( this ) result( output ) diff --git a/src/integralsTransformation/IntegralTransformation.f90 b/src/integralsTransformation/IntegralTransformation.f90 index af52afae..1e761928 100644 --- a/src/integralsTransformation/IntegralTransformation.f90 +++ b/src/integralsTransformation/IntegralTransformation.f90 @@ -173,10 +173,10 @@ program IntegralsTransformation nameOfSpecies = trim( MolecularSystem_getNameOfSpecie( i ) ) !! For PT = 2 there is no need to transform integrals for all species" - if ( partialTransform == "PT2" .and. CONTROL_instance%IONIZE_SPECIE(1) /= "NONE" ) then + if ( partialTransform == "PT2" .and. CONTROL_instance%IONIZE_SPECIES(1) /= "NONE" ) then transformThisSpecies = .False. - do z = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if ( nameOfSpecies == CONTROL_instance%IONIZE_SPECIE(z) ) then + do z = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if ( nameOfSpecies == CONTROL_instance%IONIZE_SPECIES(z) ) then transformThisSpecies = .True. end if end do @@ -256,11 +256,11 @@ program IntegralsTransformation nameOfOtherSpecie= trim( MolecularSystem_getNameOfSpecie( j ) ) !! For PT = 2 there is no need to transform integrals for all species" - if ( partialTransform == "PT2" .and. CONTROL_instance%IONIZE_SPECIE(1) /= "NONE" ) then + if ( partialTransform == "PT2" .and. CONTROL_instance%IONIZE_SPECIES(1) /= "NONE" ) then transformTheseSpecies = .False. - do z = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if ( nameOfSpecies == CONTROL_instance%IONIZE_SPECIE(z) .or. & - nameOfOtherSpecie == CONTROL_instance%IONIZE_SPECIE(z) ) then + do z = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if ( nameOfSpecies == CONTROL_instance%IONIZE_SPECIES(z) .or. & + nameOfOtherSpecie == CONTROL_instance%IONIZE_SPECIES(z) ) then transformTheseSpecies = .True. end if end do diff --git a/src/integralsTransformation/TransformIntegralsC.f90 b/src/integralsTransformation/TransformIntegralsC.f90 index 61f47ab1..ad0d62eb 100644 --- a/src/integralsTransformation/TransformIntegralsC.f90 +++ b/src/integralsTransformation/TransformIntegralsC.f90 @@ -1504,7 +1504,7 @@ subroutine TransformIntegralsC_checkMOIntegralType(speciesID, this, symmetric) if ( trim(this%partialTransform)=="PT2" ) then symmetric = .false. - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then !! all this%p_l = totalOccupation !! HOMO this%p_u = totalOccupation + 1 !! LUMO @@ -1541,8 +1541,8 @@ subroutine TransformIntegralsC_checkMOIntegralType(speciesID, this, symmetric) else if (CONTROL_instance%PT_TRANSITION_OPERATOR) then - this%p_l = CONTROL_instance%IONIZE_MO - this%p_u = CONTROL_instance%IONIZE_MO + this%p_l = CONTROL_instance%IONIZE_MO(1) + this%p_u = CONTROL_instance%IONIZE_MO(1) this%q_l = coreOrbitals+1 this%q_u = totalActiveOrbitals @@ -1552,8 +1552,8 @@ subroutine TransformIntegralsC_checkMOIntegralType(speciesID, this, symmetric) this%s_u = totalActiveOrbitals else - this%p_l = CONTROL_instance%IONIZE_MO - this%p_u = CONTROL_instance%IONIZE_MO + this%p_l = CONTROL_instance%IONIZE_MO(1) + this%p_u = CONTROL_instance%IONIZE_MO(1) this%q_l = coreOrbitals+1 this%q_u = totalActiveOrbitals @@ -1571,7 +1571,7 @@ subroutine TransformIntegralsC_checkMOIntegralType(speciesID, this, symmetric) if ( trim(this%partialTransform)=="MP2-PT2" ) then symmetric = .false. - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then this%p_l = coreOrbitals+1 this%p_u = totalOccupation + 1 @@ -1586,7 +1586,7 @@ subroutine TransformIntegralsC_checkMOIntegralType(speciesID, this, symmetric) if (CONTROL_instance%PT_TRANSITION_OPERATOR) then this%p_l = coreOrbitals+1 - this%p_u = max(CONTROL_instance%IONIZE_MO,totalOccupation) + this%p_u = max(CONTROL_instance%IONIZE_MO(1),totalOccupation) this%q_l = coreOrbitals+1 this%q_u = totalActiveOrbitals @@ -1597,7 +1597,7 @@ subroutine TransformIntegralsC_checkMOIntegralType(speciesID, this, symmetric) else this%p_l = coreOrbitals+1 - this%p_u = max(CONTROL_instance%IONIZE_MO,totalOccupation) + this%p_u = max(CONTROL_instance%IONIZE_MO(1),totalOccupation) this%q_l = coreOrbitals+1 this%q_u = totalActiveOrbitals @@ -1710,7 +1710,7 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI this%s_upperOrbital = otherTotalActiveOrbitals symmetric = .true. - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE" ) then + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE" ) then symmetric = .false. ionizeA = .false. @@ -1719,16 +1719,16 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI nameOfSpecies= trim( MolecularSystem_getNameOfSpecie( speciesID ) ) nameOfOtherSpecies= trim( MolecularSystem_getNameOfSpecie( otherSpeciesID ) ) - do s = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + do s = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeA = .true. end if - if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeB = .true. end if end do - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then if ( ionizeA .and. ionizeB ) then @@ -1766,10 +1766,10 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI end if - else if ( CONTROL_instance%IONIZE_MO /= 0 ) then !!occ and vir.. + else if ( CONTROL_instance%IONIZE_MO(1) /= 0 ) then !!occ and vir.. if ( ionizeA .and. ionizeB ) then - if ( CONTROL_instance%IONIZE_MO <= totalOccupation .and. CONTROL_instance%IONIZE_MO <= othertotalOccupation ) then + if ( CONTROL_instance%IONIZE_MO(1) <= totalOccupation .and. CONTROL_instance%IONIZE_MO(1) <= othertotalOccupation ) then this%p_lowerOrbital = coreOrbitals+1 this%p_upperOrbital = totalOccupation this%q_lowerOrbital = coreOrbitals+1 @@ -1779,14 +1779,14 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI this%s_lowerOrbital = otherCoreOrbitals+1 this%s_upperOrbital = otherTotalActiveOrbitals - else if ( CONTROL_instance%IONIZE_MO > totalOccupation .and. CONTROL_instance%IONIZE_MO > othertotalOccupation ) then + else if ( CONTROL_instance%IONIZE_MO(1) > totalOccupation .and. CONTROL_instance%IONIZE_MO(1) > othertotalOccupation ) then this%p_lowerOrbital = coreOrbitals+1 - this%p_upperOrbital = CONTROL_instance%IONIZE_MO + this%p_upperOrbital = CONTROL_instance%IONIZE_MO(1) this%q_lowerOrbital = coreOrbitals+1 this%q_upperOrbital = totalActiveOrbitals this%r_lowerOrbital = otherCoreOrbitals+1 - this%r_upperOrbital = CONTROL_instance%IONIZE_MO + this%r_upperOrbital = CONTROL_instance%IONIZE_MO(1) this%s_lowerOrbital = otherCoreOrbitals+1 this%s_upperOrbital = otherTotalActiveOrbitals @@ -1794,8 +1794,8 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI else if ( ionizeA .and. .not. ionizeB ) then - this%p_lowerOrbital = CONTROL_instance%IONIZE_MO - this%p_upperOrbital = CONTROL_instance%IONIZE_MO + this%p_lowerOrbital = CONTROL_instance%IONIZE_MO(1) + this%p_upperOrbital = CONTROL_instance%IONIZE_MO(1) this%q_lowerOrbital = coreOrbitals+1 this%q_upperOrbital = totalActiveOrbitals @@ -1810,8 +1810,8 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI this%p_upperOrbital = totalOccupation this%q_lowerOrbital = totalOccupation + 1 this%q_upperOrbital = totalActiveOrbitals - this%r_lowerOrbital = CONTROL_instance%IONIZE_MO - this%r_upperOrbital = CONTROL_instance%IONIZE_MO + this%r_lowerOrbital = CONTROL_instance%IONIZE_MO(1) + this%r_upperOrbital = CONTROL_instance%IONIZE_MO(1) this%s_lowerOrbital = otherCoreOrbitals+1 this%s_upperOrbital = otherTotalActiveOrbitals @@ -1837,7 +1837,7 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI this%s_upperOrbital = otherTotalActiveOrbitals symmetric = .true. - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE" ) then + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE" ) then symmetric = .false. ionizeA = .false. @@ -1846,16 +1846,16 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI nameOfSpecies= trim( MolecularSystem_getNameOfSpecie( speciesID ) ) nameOfOtherSpecies= trim( MolecularSystem_getNameOfSpecie( otherSpeciesID ) ) - do s = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + do s = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeA = .true. end if - if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeB = .true. end if end do - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then if ( ionizeA .and. ionizeB ) then @@ -1893,10 +1893,10 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI end if - else if ( CONTROL_instance%IONIZE_MO /= 0 ) then !!occ and vir.. + else if ( CONTROL_instance%IONIZE_MO(1) /= 0 ) then !!occ and vir.. if ( ionizeA .and. ionizeB ) then - if ( CONTROL_instance%IONIZE_MO <= totalOccupation .and. CONTROL_instance%IONIZE_MO <= othertotalOccupation ) then + if ( CONTROL_instance%IONIZE_MO(1) <= totalOccupation .and. CONTROL_instance%IONIZE_MO(1) <= othertotalOccupation ) then this%p_lowerOrbital = coreOrbitals+1 this%p_upperOrbital = totalOccupation this%q_lowerOrbital = coreOrbitals+1 @@ -1906,14 +1906,14 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI this%s_lowerOrbital = otherCoreOrbitals+1 this%s_upperOrbital = otherTotalActiveOrbitals - else if ( CONTROL_instance%IONIZE_MO > totalOccupation .and. CONTROL_instance%IONIZE_MO > othertotalOccupation ) then + else if ( CONTROL_instance%IONIZE_MO(1) > totalOccupation .and. CONTROL_instance%IONIZE_MO(1) > othertotalOccupation ) then this%p_lowerOrbital = coreOrbitals+1 - this%p_upperOrbital = CONTROL_instance%IONIZE_MO + this%p_upperOrbital = CONTROL_instance%IONIZE_MO(1) this%q_lowerOrbital = coreOrbitals+1 this%q_upperOrbital = totalActiveOrbitals this%r_lowerOrbital = otherCoreOrbitals+1 - this%r_upperOrbital = CONTROL_instance%IONIZE_MO + this%r_upperOrbital = CONTROL_instance%IONIZE_MO(1) this%s_lowerOrbital = otherCoreOrbitals+1 this%s_upperOrbital = otherTotalActiveOrbitals @@ -1922,7 +1922,7 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI else if ( ionizeA .and. .not. ionizeB ) then this%p_lowerOrbital = coreOrbitals+1 - this%p_upperOrbital = max(CONTROL_instance%IONIZE_MO,totalOccupation) + this%p_upperOrbital = max(CONTROL_instance%IONIZE_MO(1),totalOccupation) this%q_lowerOrbital = coreOrbitals+1 this%q_upperOrbital = totalActiveOrbitals @@ -1938,7 +1938,7 @@ subroutine TransformIntegralsC_checkInterMOIntegralType(speciesID, otherSpeciesI this%q_lowerOrbital = totalOccupation + 1 this%q_upperOrbital = totalActiveOrbitals this%r_lowerOrbital = otherCoreOrbitals+1 - this%r_upperOrbital = max(CONTROL_instance%IONIZE_MO,otherTotalOccupation) + this%r_upperOrbital = max(CONTROL_instance%IONIZE_MO(1),otherTotalOccupation) this%s_lowerOrbital = otherCoreOrbitals+1 this%s_upperOrbital = otherTotalActiveOrbitals diff --git a/src/integralsTransformation/TransformIntegralsE.f90 b/src/integralsTransformation/TransformIntegralsE.f90 index a390d5e6..e001e01e 100644 --- a/src/integralsTransformation/TransformIntegralsE.f90 +++ b/src/integralsTransformation/TransformIntegralsE.f90 @@ -1951,7 +1951,7 @@ subroutine TransformIntegralsE_checkMOIntegralType(speciesID, this) !! only the (ip|aq) integrals will be transformed if ( trim(this%partialTransform)=="PT2" ) then - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then !! all this%q_l = coreOrbitals+1!totalOccupation !! HOMO this%q_u = totalOccupation+1 !! LUMO @@ -1988,8 +1988,8 @@ subroutine TransformIntegralsE_checkMOIntegralType(speciesID, this) else if (CONTROL_instance%PT_TRANSITION_OPERATOR) then - this%q_l = CONTROL_instance%IONIZE_MO - this%q_u = CONTROL_instance%IONIZE_MO + this%q_l = CONTROL_instance%IONIZE_MO(1) + this%q_u = CONTROL_instance%IONIZE_MO(1) this%p_l = coreOrbitals+1 this%p_u = totalActiveOrbitals @@ -2001,8 +2001,8 @@ subroutine TransformIntegralsE_checkMOIntegralType(speciesID, this) else - this%q_l = CONTROL_instance%IONIZE_MO - this%q_u = CONTROL_instance%IONIZE_MO + this%q_l = CONTROL_instance%IONIZE_MO(1) + this%q_u = CONTROL_instance%IONIZE_MO(1) this%p_l = coreOrbitals+1 this%p_u = totalActiveOrbitals @@ -2019,7 +2019,7 @@ subroutine TransformIntegralsE_checkMOIntegralType(speciesID, this) !for a simultaneous PT2 - MP2 calculation if ( trim(this%partialTransform)=="MP2-PT2" ) then - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then !! all this%q_l = coreOrbitals+1!totalOccupation !! HOMO this%q_u = totalOccupation + 1 !! LUMO @@ -2035,7 +2035,7 @@ subroutine TransformIntegralsE_checkMOIntegralType(speciesID, this) if (CONTROL_instance%PT_TRANSITION_OPERATOR) then this%q_l = coreOrbitals+1 - this%q_u = max(CONTROL_instance%IONIZE_MO,totalOccupation) + this%q_u = max(CONTROL_instance%IONIZE_MO(1),totalOccupation) this%p_l = coreOrbitals+1 this%p_u = totalActiveOrbitals @@ -2048,7 +2048,7 @@ subroutine TransformIntegralsE_checkMOIntegralType(speciesID, this) else this%q_l = coreOrbitals+1 - this%q_u = max(CONTROL_instance%IONIZE_MO,totalOccupation) + this%q_u = max(CONTROL_instance%IONIZE_MO(1),totalOccupation) this%p_l = coreOrbitals+1 this%p_u = totalActiveOrbitals @@ -2144,7 +2144,7 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI this%r_l = coreOrbitals+1 this%r_u = otherTotalActiveOrbitals - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE" ) then + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE" ) then ionizeA = .false. ionizeB = .false. @@ -2152,16 +2152,16 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI nameOfSpecies= trim( MolecularSystem_getNameOfSpecie( speciesID ) ) nameOfOtherSpecies= trim( MolecularSystem_getNameOfSpecie( otherSpeciesID ) ) - do s = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + do s = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeA = .true. end if - if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeB = .true. end if end do - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then if ( ionizeA .and. ionizeB ) then @@ -2199,10 +2199,10 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI end if - else if ( CONTROL_instance%IONIZE_MO /= 0 ) then !!occ and vir.. + else if ( CONTROL_instance%IONIZE_MO(1) /= 0 ) then !!occ and vir.. if ( ionizeA .and. ionizeB ) then - if ( CONTROL_instance%IONIZE_MO <= totalOccupation .and. CONTROL_instance%IONIZE_MO <= othertotalOccupation ) then + if ( CONTROL_instance%IONIZE_MO(1) <= totalOccupation .and. CONTROL_instance%IONIZE_MO(1) <= othertotalOccupation ) then this%q_l = coreOrbitals+1 this%q_u = totalOccupation this%p_l = coreOrbitals+1 @@ -2212,14 +2212,14 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI this%r_l = otherCoreOrbitals+1 this%r_u = otherTotalActiveOrbitals - else if ( CONTROL_instance%IONIZE_MO > totalOccupation .and. CONTROL_instance%IONIZE_MO > othertotalOccupation ) then + else if ( CONTROL_instance%IONIZE_MO(1) > totalOccupation .and. CONTROL_instance%IONIZE_MO(1) > othertotalOccupation ) then this%q_l = coreOrbitals+1 - this%q_u = totalActiveOrbitals!CONTROL_instance%IONIZE_MO + this%q_u = totalActiveOrbitals!CONTROL_instance%IONIZE_MO(1) this%p_l = coreOrbitals+1 this%p_u = totalActiveOrbitals this%s_l = otherCoreOrbitals+1 - this%s_u = otherTotalActiveOrbitals!CONTROL_instance%IONIZE_MO + this%s_u = otherTotalActiveOrbitals!CONTROL_instance%IONIZE_MO(1) this%r_l = otherCoreOrbitals+1 this%r_u = otherTotalActiveOrbitals @@ -2227,8 +2227,8 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI else if ( ionizeA .and. .not. ionizeB ) then - this%q_l = CONTROL_instance%IONIZE_MO - this%q_u = CONTROL_instance%IONIZE_MO + this%q_l = CONTROL_instance%IONIZE_MO(1) + this%q_u = CONTROL_instance%IONIZE_MO(1) if (CONTROL_instance%PT_TRANSITION_OPERATOR) then @@ -2250,8 +2250,8 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI this%p_l = totalOccupation + 1 this%p_u = totalActiveOrbitals - !this%s_l = CONTROL_instance%IONIZE_MO !... - !this%s_u = CONTROL_instance%IONIZE_MO !... + !this%s_l = CONTROL_instance%IONIZE_MO(1) !... + !this%s_u = CONTROL_instance%IONIZE_MO(1) !... this%s_l = otherCoreOrbitals+1 this%s_u = otherTotalActiveOrbitals @@ -2264,7 +2264,7 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI end if - !if ( CONTROL_instance%IONIZE_MO /= 0 ) then + !if ( CONTROL_instance%IONIZE_MO(1) /= 0 ) then end if end if @@ -2282,7 +2282,7 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI this%r_l = otherCoreOrbitals+1 this%r_u = totalActiveOrbitals - if (CONTROL_instance%IONIZE_SPECIE(1) /= "NONE" ) then + if (CONTROL_instance%IONIZE_SPECIES(1) /= "NONE" ) then ionizeA = .false. ionizeB = .false. @@ -2290,16 +2290,16 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI nameOfSpecies= trim( MolecularSystem_getNameOfSpecie( speciesID ) ) nameOfOtherSpecies= trim( MolecularSystem_getNameOfSpecie( otherSpeciesID ) ) - do s = 1, size(CONTROL_instance%IONIZE_SPECIE ) - if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + do s = 1, size(CONTROL_instance%IONIZE_SPECIES ) + if ( nameOfSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeA = .true. end if - if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIE(s)) ) then + if ( nameOfOtherSpecies == trim(CONTROL_instance%IONIZE_SPECIES(s)) ) then ionizeB = .true. end if end do - if ( CONTROL_instance%IONIZE_MO == 0 ) then + if ( CONTROL_instance%IONIZE_MO(1) == 0 ) then if ( ionizeA .and. ionizeB ) then @@ -2337,10 +2337,10 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI end if - else if ( CONTROL_instance%IONIZE_MO /= 0 ) then !!occ and vir.. + else if ( CONTROL_instance%IONIZE_MO(1) /= 0 ) then !!occ and vir.. if ( ionizeA .and. ionizeB ) then - if ( CONTROL_instance%IONIZE_MO <= totalOccupation .and. CONTROL_instance%IONIZE_MO <= othertotalOccupation ) then + if ( CONTROL_instance%IONIZE_MO(1) <= totalOccupation .and. CONTROL_instance%IONIZE_MO(1) <= othertotalOccupation ) then this%q_l = coreOrbitals+1 this%q_u = totalOccupation this%p_l = coreOrbitals+1 @@ -2350,14 +2350,14 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI this%r_l = otherCoreOrbitals+1 this%r_u = otherTotalActiveOrbitals - else if ( CONTROL_instance%IONIZE_MO > totalOccupation .and. CONTROL_instance%IONIZE_MO > othertotalOccupation ) then + else if ( CONTROL_instance%IONIZE_MO(1) > totalOccupation .and. CONTROL_instance%IONIZE_MO(1) > othertotalOccupation ) then this%q_l = coreOrbitals+1 - this%q_u = totalActiveOrbitals!CONTROL_instance%IONIZE_MO + this%q_u = totalActiveOrbitals!CONTROL_instance%IONIZE_MO(1) this%p_l = coreOrbitals+1 this%p_u = totalActiveOrbitals this%s_l = otherCoreOrbitals+1 - this%s_u = otherTotalActiveOrbitals!CONTROL_instance%IONIZE_MO + this%s_u = otherTotalActiveOrbitals!CONTROL_instance%IONIZE_MO(1) this%r_l = otherCoreOrbitals+1 this%r_u = otherTotalActiveOrbitals @@ -2366,7 +2366,7 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI else if ( ionizeA .and. .not. ionizeB ) then this%q_l = coreOrbitals+1 - this%q_u = max(CONTROL_instance%IONIZE_MO,totalOccupation) + this%q_u = max(CONTROL_instance%IONIZE_MO(1),totalOccupation) if (CONTROL_instance%PT_TRANSITION_OPERATOR) then this%q_l = coreOrbitals+1 @@ -2387,8 +2387,8 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI this%p_l = totalOccupation + 1 this%p_u = totalActiveOrbitals - !this%s_l = CONTROL_instance%IONIZE_MO !... - !this%s_u = CONTROL_instance%IONIZE_MO !... + !this%s_l = CONTROL_instance%IONIZE_MO(1) !... + !this%s_u = CONTROL_instance%IONIZE_MO(1) !... this%s_l = otherCoreOrbitals+1 this%s_u = otherTotalActiveOrbitals @@ -2401,7 +2401,7 @@ subroutine TransformIntegralsE_checkInterMOIntegralType(speciesID, otherSpeciesI end if - !if ( CONTROL_instance%IONIZE_MO /= 0 ) then + !if ( CONTROL_instance%IONIZE_MO(1) /= 0 ) then end if diff --git a/src/ints/Ints.f90 b/src/ints/Ints.f90 index 740d2fb9..9090bbaf 100644 --- a/src/ints/Ints.f90 +++ b/src/ints/Ints.f90 @@ -57,9 +57,6 @@ Program Ints !!Load CONTROL Parameters call MolecularSystem_loadFromFile( "LOWDIN.DAT" ) - !!Load the system in lowdin.bas format - call MolecularSystem_loadFromFile( "LOWDIN.BAS" ) - !!Load the system in lowdin.sys format call MolecularSystem_loadFromFile( "LOWDIN.SYS" ) diff --git a/src/output/InputOutput.f90 b/src/output/InputOutput.f90 index 73ad9f6d..80513061 100644 --- a/src/output/InputOutput.f90 +++ b/src/output/InputOutput.f90 @@ -169,11 +169,13 @@ subroutine InputOutput_load( ) Output_type="" Output_species="ALL" Output_state=1 - Output_orbital=0 - Output_dimensions=0 - Output_cubeSize=0.0_8 + Output_orbital=1 + Output_dimensions=2 + Output_cubeSize=5.0_8 Output_point1(:)=0.0_8 + Output_point1(3)=-5.0_8 Output_point2(:)=0.0_8 + Output_point2(3)=5.0_8 Output_point3(:)=0.0_8 read(4,NML=Output, iostat=stat) diff --git a/src/output/OutputBuilder.f90 b/src/output/OutputBuilder.f90 index 5e13fdda..8637f48b 100644 --- a/src/output/OutputBuilder.f90 +++ b/src/output/OutputBuilder.f90 @@ -125,32 +125,40 @@ subroutine OutputBuilder_constructor(this, ID, type ,species, state, orbital, di ! print *, "this%type", this%type this%outputID=ID this%species=trim(String_getUppercase(species)) - if( trim(this%species) .eq. "ALL" ) then + if( trim(this%species) .eq. "ALL" .and. this%type .ne. "ORBITALPLOT" ) then allocate(this%fileName(MolecularSystem_getNumberOfQuantumSpecies())) + else if( trim(this%species) .eq. "ALL" .and. this%type .eq. "ORBITALPLOT" ) then + allocate(this%fileName(1)) + this%species=MolecularSystem_getNameOfSpecie(1) else allocate(this%fileName(1)) end if + this%state=1 if( present(state)) this%state=state - this%orbital=0 + this%orbital=1 if( present(orbital)) this%orbital=orbital - this%dimensions=0 + this%dimensions=2 if( present(dimensions)) this%dimensions=dimensions - this%cubeSize=0 + this%cubeSize=10 if( present(cubeSize)) this%cubeSize=cubeSize call Vector_constructor(this%point1, 3, 0.0_8 ) call Vector_constructor(this%point2, 3, 0.0_8 ) call Vector_constructor(this%point3, 3, 0.0_8 ) + this%point1%values(3)=-5.0 + this%point2%values(3)=5.0 + if( present(point1)) this%point1%values=point1%values if( present(point2)) this%point2%values=point2%values if( present(point3)) this%point3%values=point3%values if ( trim(CONTROL_instance%UNITS) == "ANGS") then - this%point1%values= this%point1%values / AMSTRONG - this%point2%values= this%point2%values / AMSTRONG - this%point3%values= this%point3%values / AMSTRONG + this%point1%values= this%point1%values/AMSTRONG + this%point2%values= this%point2%values/AMSTRONG + this%point3%values= this%point3%values/AMSTRONG + this%cubeSize=this%cubeSize/AMSTRONG end if this%auxID=1 @@ -228,15 +236,49 @@ subroutine OutputBuilder_show(this) end do ! TODO Fix this line. ! if (this%filename2 /= "") print *, "FileName 2: ", this%fileName2 - if (this%species /= "ALL") write (*,"(A20,A10)") "for species: ", this%species - if (this%orbital /= 0) write (*,"(A20,I10)") "for orbital: ", this%orbital + if (this%species /= "ALL") write (*,"(A20,A)") "for species: ", trim(this%species) if (this%state /= 1) write (*,"(A20,I10)") "for excited state: ", this%state - if (this%dimensions /= 0) write (*,"(A20,I2)") "dimensions: ", this%dimensions - if (this%cubeSize /= 0.0_8) write (*,"(A20,F15.5)") "cube size (a.u.): ", this%cubeSize - if (this%cubeSize /= 0.0_8) write (*,"(A20,F15.5)") "cube center (a.u.): ", this%point1%values(1) - if (this%dimensions >= 1) write (*,"(A20,F10.5,F10.5,F10.5)") "Point 1 (a.u.): ", this%point1%values(1), this%point1%values(2), this%point1%values(3) - if (this%dimensions >= 2) write (*,"(A20,F10.5,F10.5,F10.5)") "Point 2 (a.u.): ", this%point2%values(1), this%point2%values(2), this%point2%values(3) - if (this%dimensions >= 3) write (*,"(A20,F10.5,F10.5,F10.5)") "Point 3 (a.u.): ", this%point3%values(1), this%point3%values(2), this%point3%values(3) + + select case(trim(this%type)) + + case ( "MOLDENFILE") + + case ("VECGAMESSFILE") + + case ("CASINOFILE") + + case ("EIGENGAMESSFILE") + + case ("FCHKFILE") + + case ( "WFNFILE") + + case ( "NBO47FILE") + + case ( "WFXFILE" ) + + case ( "EXTENDEDWFNFILE") + + case ( "DENSITYPLOT") + write (*,"(A20,I10)") "dimensions: ", this%dimensions + write (*,"(A20,F10.5,F10.5,F10.5)") "Point 1 (a.u.): ", this%point1%values(1), this%point1%values(2), this%point1%values(3) + write (*,"(A20,F10.5,F10.5,F10.5)") "Point 2 (a.u.): ", this%point2%values(1), this%point2%values(2), this%point2%values(3) + if (this%dimensions >= 3) write (*,"(A20,F10.5,F10.5,F10.5)") "Point 3 (a.u.): ", this%point3%values(1), this%point3%values(2), this%point3%values(3) + + case ( "DENSITYCUBE") + write (*,"(A20,F10.5)") "cube size (a.u.): ", this%cubeSize + write (*,"(A20,3F10.5)") "cube center (a.u.): ", this%point1%values(1:3) + + case ( "ORBITALPLOT") + write (*,"(A20,I10)") "for orbital: ", this%orbital + write (*,"(A20,I10)") "dimensions: ", this%dimensions + write (*,"(A20,F10.5,F10.5,F10.5)") "Point 1 (a.u.): ", this%point1%values(1), this%point1%values(2), this%point1%values(3) + write (*,"(A20,F10.5,F10.5,F10.5)") "Point 2 (a.u.): ", this%point2%values(1), this%point2%values(2), this%point2%values(3) + if (this%dimensions >= 3) write (*,"(A20,F10.5,F10.5,F10.5)") "Point 3 (a.u.): ", this%point3%values(1), this%point3%values(2), this%point3%values(3) + + case default + + end select print *, "--------------------------------------------------------" print *, "" @@ -1812,11 +1854,11 @@ subroutine OutputBuilder_get3DPlot(this) if (this%type .eq. "FUKUIPLOT") write(11,*) "" do j=0,numberOfSteps n=n+1 - write (10,"(T10,F20.8,F20.8,F20.8)") i*Vector_norm(step1),j*Vector_norm(step2),val%values(n) + write (10,"(T10,F20.8,F20.8,ES20.8)") i*Vector_norm(step1),j*Vector_norm(step2),val%values(n) if (val%values(n) > maxValue) maxValue = val%values(n) if (val%values(n) < minValue) minValue = val%values(n) if (this%type .eq. "FUKUIPLOT" ) then - write (11,"(T10,F20.8,F20.8,F20.8)") i*Vector_norm(step1),j*Vector_norm(step2),val2%values(n) + write (11,"(T10,F20.8,F20.8,ES20.8)") i*Vector_norm(step1),j*Vector_norm(step2),val2%values(n) if (val2%values(n) > maxValue2) maxValue2 = val2%values(n) if (val2%values(n) < minValue2) minValue2 = val2%values(n) end if @@ -1907,8 +1949,8 @@ subroutine OutputBuilder_get2DPlot(this) n=0 do i=0,numberOfSteps n=n+1 - write (10,"(T10,F20.8,F20.8)") i*Vector_norm(step),val%values(n) - if (this%type .eq. "FUKUIPLOT") write (11,"(T10,F20.8,F20.8)") i*Vector_norm(step),val2%values(n) + write (10,"(T10,F20.8,ES20.8)") i*Vector_norm(step),val%values(n) + if (this%type .eq. "FUKUIPLOT") write (11,"(T10,F20.8,ES20.8)") i*Vector_norm(step),val2%values(n) end do close(10) @@ -2166,8 +2208,13 @@ subroutine OutputBuilder_getDensityPlot(this) write (10,*) "" do j=0,numberOfSteps n=n+1 - write (10,"(T10,F20.8,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1), -plotDistance2*0.5+j*Vector_norm(step2), & - val%values(n) + if(val%values(n)>1.0E-99_8) then + write (10,"(T10,F20.8,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1), -plotDistance2*0.5+j*Vector_norm(step2), & + val%values(n) + else + write (10,"(T10,F20.8,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1), -plotDistance2*0.5+j*Vector_norm(step2), & + 0.0 + end if if (val%values(n) > maxValue) maxValue = val%values(n) if (val%values(n) < minValue) minValue = val%values(n) ! print *, coordinate, val @@ -2204,7 +2251,11 @@ subroutine OutputBuilder_getDensityPlot(this) n=0 do i=0,numberOfSteps n=n+1 - write (10,"(T10,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1),val%values(n) + if(val%values(n)>1.0E-99_8) then + write (10,"(T10,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1),val%values(n) + else + write (10,"(T10,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1),0.0 + end if ! print *, coordinate, val end do @@ -2269,14 +2320,14 @@ subroutine OutputBuilder_make2DGraph(fileName, title, x_title, y_title,& write (10,"(A)") 'set ylabel "'//trim(y_title)//'"' write (10,"(A)") 'set format y "'//trim(auxYformat)//'"' if( auxNumOfGraphs >1) then - write (10,"(A$)") 'plot '//trim(auxXRange)//trim(auxYRange)//' "'//trim(fileName)//'" using 1:2 w l title "" smooth csplines' + write (10,"(A$)") 'plot '//trim(auxXRange)//trim(auxYRange)//' "'//trim(fileName)//'" using 1:2 w l title "" ' do i=2, auxNumOfGraphs charNumOfGraph=String_convertIntegerToString(i+1) - write (10,"(A$)") ', "'//trim(fileName)//'.dat"'//' using 1:'//trim(charNumOfGraph)//' w l title "" smooth csplines' + write (10,"(A$)") ', "'//trim(fileName)//'.dat"'//' using 1:'//trim(charNumOfGraph)//' w l title "" ' end do write (10,"(A)") "" else - write (10,"(A)") 'plot '//trim(auxXRange)//trim(auxYRange)//' "'//trim(fileName)//'" w l title "" smooth csplines' + write (10,"(A)") 'plot '//trim(auxXRange)//trim(auxYRange)//' "'//trim(fileName)//'" w l title "" ' end if write (10,"(A)") 'set output' close(10) diff --git a/src/scf/DensityMatrixSCFGuess.f90 b/src/scf/DensityMatrixSCFGuess.f90 index 7b274eae..a4937b0a 100644 --- a/src/scf/DensityMatrixSCFGuess.f90 +++ b/src/scf/DensityMatrixSCFGuess.f90 @@ -50,6 +50,7 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio type(Matrix), intent(inout) :: orbitals logical, intent(in) :: printInfo + type(Matrix) :: auxMatrix character(30) :: nameOfSpecies integer(8) :: orderOfMatrix, occupationNumber logical :: existPlain, existBinnary, readSuccess @@ -132,23 +133,38 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio call DensityMatrixSCFGuess_exception( ERROR, "the selected guess method for "//nameOfSpecies//" is not implemented", "at program SCF module DensityMatrixSCFGuess") end select - end if + if(CONTROL_instance%DEBUG_SCFS) then + print *, "Guess orbitals for", nameOfSpecies + call Matrix_show(orbitals) + end if + + call Matrix_copyConstructor(auxMatrix,orbitals) + !! Segment for fractional occupations: introduce fractional occupation + if (trim(nameOfSpecies) == trim(CONTROL_instance%IONIZE_SPECIES(1)) ) then + do i=1,size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .lt. 1.0_8) then + if(printInfo) write (*,"(A,F6.2,A,I5,A,A)") "Removing ", (1.0-CONTROL_instance%MO_FRACTION_OCCUPATION(i))*100, & + " % of the density associated with orbital No. ", CONTROL_instance%IONIZE_MO(i), " of ", trim(nameOfSpecies) + auxMatrix%values(:,CONTROL_instance%IONIZE_MO(i)) = auxMatrix%values(:,CONTROL_instance%IONIZE_MO(i))*sqrt(CONTROL_instance%MO_FRACTION_OCCUPATION(i)) + end if + end do + end if + do i = 1 , orderOfMatrix do j = 1 , orderOfMatrix do k = 1 , occupationNumber - densityMatrix%values(i,j) = densityMatrix%values( i,j ) + orbitals%values(i,k) * orbitals%values(j,k) + densityMatrix%values(i,j) = densityMatrix%values( i,j ) + auxMatrix%values(i,k) * auxMatrix%values(j,k) end do end do end do densityMatrix%values=densityMatrix%values*MolecularSystem_getEta( speciesID ) - + if ( CONTROL_instance%BUILD_MIXED_DENSITY_MATRIX ) then densityMatrix%values(occupationNumber,:) = 0.1*densityMatrix%values(occupationNumber,:)*densityMatrix%values(occupationNumber+1,:) end if - end subroutine DensityMatrixSCFGuess_getGuess !> diff --git a/src/scf/MultiSCF.f90 b/src/scf/MultiSCF.f90 index 2b502608..4cc3be71 100644 --- a/src/scf/MultiSCF.f90 +++ b/src/scf/MultiSCF.f90 @@ -618,8 +618,8 @@ subroutine MultiSCF_getInitialGuess(this,wfObjects) type(MultiSCF) :: this type(WaveFunction) :: wfObjects(*) - integer :: speciesID, otherSpeciesID - real(8) :: normCheck + integer :: speciesID, otherSpeciesID, i + real(8) :: expectedOccupation,normCheck !!********************************************************** !! Build Guess and first density matrix @@ -636,12 +636,20 @@ subroutine MultiSCF_getInitialGuess(this,wfObjects) write(*,"(A15,A10,A40,F12.6)") "number of ", trim(MolecularSystem_getNameOfSpecies( speciesID )) , & " particles in guess density matrix: ", normCheck - if(abs(MolecularSystem_getEta(speciesID)*MolecularSystem_instance%species(speciesID)%ocupationNumber/normCheck-1.0) .gt. 1.0E-2 ) then - wfObjects(speciesID)%densityMatrix%values=wfObjects(speciesID)%densityMatrix%values*& - MolecularSystem_getEta(speciesID)*MolecularSystem_instance%species(speciesID)%ocupationNumber/normCheck + expectedOccupation=MolecularSystem_getEta(speciesID)*MolecularSystem_instance%species(speciesID)%ocupationNumber + if (trim(MolecularSystem_getNameOfSpecies( speciesID )) .eq. trim(CONTROL_instance%IONIZE_SPECIES(1))) then + do i=1,size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .lt. 1.0_8) & + expectedOccupation=expectedOccupation-MolecularSystem_getEta(speciesID)*(1.0-CONTROL_instance%MO_FRACTION_OCCUPATION(i)) + end do + end if + + if (abs(expectedOccupation/normCheck-1.0) .gt. 1.0E-3 ) then + ! wfObjects(speciesID)%densityMatrix%values=wfObjects(speciesID)%densityMatrix%values*expectedOccupation/normCheck if ( this%printSCFiterations ) & - write(*,"(A65,F12.6)") "renormalized density matrix with deviation from number of particles", & - MolecularSystem_getEta(speciesID)*MolecularSystem_instance%species(speciesID)%ocupationNumber + write(*,"(A65,F12.6)") "Warning!, density matrix with deviation from number of particles", expectedOccupation + ! write(*,"(A65,F12.6)") "renormalized density matrix with deviation from number of particles", & + ! sum( transpose(wfObjects(speciesID)%densityMatrix%values)*wfObjects(speciesID)%overlapMatrix%values) end if if ( CONTROL_instance%DEBUG_SCFS ) then @@ -799,7 +807,6 @@ subroutine MultiSCF_solveHartreeFockRoothan(this,wfObjects,libint2Objects) ! print *,"" ! print *,"...end Second Multi-Species SCF calculation" ! print *,"" - ! end if call MultiSCF_obtainFinalEnergy(this,wfObjects,libint2Objects) @@ -815,7 +822,7 @@ subroutine MultiSCF_obtainFinalEnergy(this,wfObjects,libint2Objects) integer :: numberOfSpecies integer :: wfnUnit, densUnit - integer :: speciesID, otherSpeciesID + integer :: speciesID, otherSpeciesID, i character(50) :: wfnFile, densFile character(30) :: labels(2) character(50) :: integralsFile @@ -873,6 +880,14 @@ subroutine MultiSCF_obtainFinalEnergy(this,wfObjects,libint2Objects) call WaveFunction_buildDensityMatrix(wfObjects(speciesID)) end do + ! Reorder coefficients after creating the last density matrix + do i=1,size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .eq. 0.0_8) then + call MultiSCF_reorderIonizedCoefficients(this,wfObjects) + exit + end if + end do + !Forces equal coefficients for E-ALPHA and E-BETA in open shell calculations if ( CONTROL_instance%FORCE_CLOSED_SHELL .and. & (CONTROL_instance%METHOD .eq. "UKS" .or. CONTROL_instance%METHOD .eq. "UHF") ) then @@ -1305,6 +1320,16 @@ subroutine MultiSCF_saveWfn(this,wfObjects) labels(1) = "DENSITY" call Matrix_writeToFile(wfObjects(speciesID)%densityMatrix, unit=wfnUnit, binary=.true., arguments = labels ) + labels(1) = "KINETIC" + call Matrix_writeToFile(wfObjects(speciesID)%kineticMatrix, unit=wfnUnit, binary=.true., arguments = labels ) + + labels(1) = "ATTRACTION" + call Matrix_writeToFile(wfObjects(speciesID)%puntualInteractionMatrix, unit=wfnUnit, binary=.true., arguments = labels ) + + labels(1) = "EXTERNAL" + if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) & + call Matrix_writeToFile(wfObjects(speciesID)%externalPotentialMatrix, unit=wfnUnit, binary=.true., arguments = labels ) + labels(1) = "HCORE" call Matrix_writeToFile(wfObjects(speciesID)%hcoreMatrix, unit=wfnUnit, binary=.true., arguments = labels ) @@ -1390,5 +1415,43 @@ subroutine MultiSCF_saveWfn(this,wfObjects) close(wfnUnit) end subroutine MultiSCF_saveWfn + +!!Reorders full ionized orbitals + subroutine MultiSCF_reorderIonizedCoefficients(this,wfObjects) + implicit none + type(MultiSCF) :: this + type(WaveFunction) :: wfObjects(*) + type(Matrix) :: auxMatrix + type(Vector) :: auxVector + integer :: occupationNumber, newOccupationNumber, i, j, speciesID + + do speciesID=1, MolecularSystem_getNumberOfQuantumSpecies() + if (trim(wfObjects(speciesID)%name) .eq. trim(CONTROL_instance%IONIZE_SPECIES(1)) ) then + occupationNumber=MolecularSystem_getOcupationNumber(speciesID) + newOccupationNumber=occupationNumber + call Matrix_copyConstructor(auxMatrix,wfObjects(speciesID)%waveFunctionCoefficients) + call Vector_copyConstructor(auxVector,wfObjects(speciesID)%molecularOrbitalsEnergy) + do i= 1, size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .eq. 0.0_8 ) then + wfObjects(speciesID)%waveFunctionCoefficients%values(:,occupationNumber)=auxMatrix%values(:,CONTROL_instance%IONIZE_MO(i)) + wfObjects(speciesID)%molecularOrbitalsEnergy%values(occupationNumber)=auxVector%values(CONTROL_instance%IONIZE_MO(i)) + do j=CONTROL_instance%IONIZE_MO(i),occupationNumber-1 + wfObjects(speciesID)%waveFunctionCoefficients%values(:,j)=auxMatrix%values(:,j+1) + wfObjects(speciesID)%molecularOrbitalsEnergy%values(j)=auxVector%values(j+1) + end do + newOccupationNumber=newOccupationNumber-1 + end if + end do + molecularSystem_instance%species(speciesID)%ocupationNumber=newOccupationNumber + if(CONTROL_instance%DEBUG_SCFS) then + print *, "newOccupationNumber for", trim(wfObjects(speciesID)%name), molecularSystem_instance%species(speciesID)%ocupationNumber + call Matrix_show(wfObjects(speciesID)%waveFunctionCoefficients) + end if + end if + end do + + end subroutine MultiSCF_reorderIonizedCoefficients + + end module MultiSCF_ diff --git a/src/scf/SingleSCF.f90 b/src/scf/SingleSCF.f90 index 084a989d..6d528083 100644 --- a/src/scf/SingleSCF.f90 +++ b/src/scf/SingleSCF.f90 @@ -149,27 +149,12 @@ subroutine SingleSCF_iterate(wfObject) type(Matrix) :: fockMatrixTransformed type(Matrix) :: auxiliaryMatrix - type(Matrix) :: previousWavefunctionCoefficients, matchingMatrix, matchingMatrix2 - type(Vector) :: deltaVector - type(Matrix) :: auxOverlapMatrix + type(Matrix) :: previousWavefunctionCoefficients - integer(8) :: total3 integer(8) :: numberOfContractions - real(8) :: threshold, hold - real(8) :: levelShiftingFactor, normCheck - integer :: ocupados, occupatedOrbitals, prodigals - integer :: totales - integer :: search, trial - integer :: ii, jj, astrayOrbitals, startIteration - integer :: i,j, mu, nu, index - logical :: existFile - - character(50) :: wfnFile - character(50) :: arguments(20) - integer :: wfnUnit + real(8) :: hold ! wfnFile = trim(CONTROL_instance%INPUT_FILE)//"lowdin.vec" - wfnUnit = 30 numberOfContractions = MolecularSystem_getTotalnumberOfContractions(wfObject%species ) !!********************************************************************************************** @@ -200,7 +185,7 @@ subroutine SingleSCF_iterate(wfObject) call Matrix_copyConstructor( fockMatrixTransformed, wfObject%fockMatrix ) - if (CONTROL_instance%MO_FRACTION_OCCUPATION < 1.0_8 .or. CONTROL_instance%EXCHANGE_ORBITALS_IN_SCF ) then + if (CONTROL_instance%IONIZE_MO(1) > 0 .or. CONTROL_instance%EXCHANGE_ORBITALS_IN_SCF ) then call Matrix_copyConstructor( previousWavefunctionCoefficients, wfObject%waveFunctionCoefficients ) end if @@ -224,39 +209,12 @@ subroutine SingleSCF_iterate(wfObject) ! call Matrix_show(wfObject%waveFunctionCoefficients) !! Filtra los orbitales eliminados por el umbral de solapamiento - ! print *, "threshold", CONTROL_instance%OVERLAP_EIGEN_THRESHOLD - i=0 - do index = 1 , numberOfContractions - i=i+1 - normCheck=0.0 - if ( abs(wfObject%molecularOrbitalsEnergy%values(i)) .lt. CONTROL_instance%OVERLAP_EIGEN_THRESHOLD ) then - do mu = 1 , numberOfContractions - do nu = 1 , numberOfContractions - normCheck=normCheck+wfObject%waveFunctionCoefficients%values(mu,i)*& - wfObject%waveFunctionCoefficients%values(nu,i)*& - wfObject%overlapMatrix%values(mu,nu) - end do - end do - ! print *, "eigenvalue", i, wfObject%molecularOrbitalsEnergy%values(i), "normCheck", normCheck - - if ( normCheck .lt. CONTROL_instance%OVERLAP_EIGEN_THRESHOLD) then - ! Shift orbital coefficients to the end of the matrix and Make energy a very large number - do j = i , numberOfContractions-1 - wfObject%molecularOrbitalsEnergy%values(j)=wfObject%molecularOrbitalsEnergy%values(j+1) - wfObject%waveFunctionCoefficients%values(:,j) = wfObject%waveFunctionCoefficients%values(:,j+1) - end do - wfObject%molecularOrbitalsEnergy%values(numberOfContractions)=1/CONTROL_instance%OVERLAP_EIGEN_THRESHOLD - wfObject%waveFunctionCoefficients%values(:,numberOfContractions)=0.0 - i=i-1 - end if - - end if - end do + call Wavefunction_removeOrbitalsBelowEigenThreshold(wfObject) if ( CONTROL_instance%ACTIVATE_LEVEL_SHIFTING) & call Convergence_removeLevelShifting(wfObject%convergenceMethod,wfObject%molecularOrbitalsEnergy) !! - !! Interation ends + !! Iteration ends !!********************************************************************************************** if ( CONTROL_instance%DEBUG_SCFS) then @@ -264,14 +222,14 @@ subroutine SingleSCF_iterate(wfObject) call Matrix_show(wfObject%waveFunctionCoefficients) end if - if (CONTROL_instance%MO_FRACTION_OCCUPATION < 1.0_8 .or. CONTROL_instance%EXCHANGE_ORBITALS_IN_SCF ) & + if (CONTROL_instance%IONIZE_MO(1) > 0 .or. CONTROL_instance%EXCHANGE_ORBITALS_IN_SCF ) & call SingleSCF_orbitalExchange(wfObject,previousWavefunctionCoefficients) if ( CONTROL_instance%NO_SCF .or. CONTROL_instance%FREEZE_NON_ELECTRONIC_ORBITALS .or. CONTROL_instance%FREEZE_ELECTRONIC_ORBITALS) & call SingleSCF_readFrozenOrbitals(wfObject) !! Exchanging orbitals just for calculation excited states - if(wfObject%name == trim(CONTROL_instance%EXCITE_SPECIE) ) then + if(wfObject%name == trim(CONTROL_instance%EXCITE_SPECIES) ) then call Matrix_constructor (auxiliaryMatrix, numberOfContractions, numberOfContractions) call Matrix_copyConstructor( auxiliaryMatrix, wfObject%waveFunctionCoefficients ) @@ -380,40 +338,59 @@ subroutine SingleSCF_exception( typeMessage, description, debugDescription) end subroutine SingleSCF_exception !! @brief Orbital exchange for TOM calculation - !! + !! Handles two different cases of orbital exchange + !! When the user explicitly requires EXCHANGE_ORBITALS_IN_SCF to have a solution with max overlap to the guess function + !! Or when an orbital is selected for partial ionization + !! Last update jan-24 Felix subroutine SingleSCF_orbitalExchange(wfObject,previousWavefunctionCoefficients) type(WaveFunction) :: wfObject type(Matrix) :: previousWavefunctionCoefficients - type(Matrix) :: fockMatrixTransformed type(Matrix) :: auxiliaryMatrix type(Matrix) :: matchingMatrix, matchingMatrix2 type(Vector) :: deltaVector + type(IVector) :: orbitalsVector type(Matrix) :: auxOverlapMatrix - integer(8) :: total3 integer(8) :: numberOfContractions - real(8) :: threshold, hold - real(8) :: levelShiftingFactor, normCheck - integer :: ocupados, occupatedOrbitals, prodigals - integer :: totales + real(8) :: threshold, maxOverlap, hold + integer :: activeOrbitals, prodigals integer :: search, trial integer :: ii, jj, astrayOrbitals, startIteration - integer :: i,j, mu, nu, index - logical :: existFile - - character(50) :: wfnFile - character(50) :: arguments(20) - integer :: wfnUnit - - startIteration=1 - threshold=0.8_8 - + integer :: i + + startIteration=0 + threshold=CONTROL_instance%EXCHANGE_ORBITAL_THRESHOLD + + !Handles two different cases of orbital exchange + !When the user explicitly requires EXCHANGE_ORBITALS_IN_SCF to have a solution with max overlap to the guess function + !Or when an orbital is selected for partial ionization + if(CONTROL_instance%EXCHANGE_ORBITALS_IN_SCF) then + activeOrbitals = MolecularSystem_getOcupationNumber(wfObject%species) + call Vector_constructorInteger(orbitalsVector,activeOrbitals) + do i=1,activeOrbitals + orbitalsVector%values(i)=i + end do + else if(CONTROL_instance%IONIZE_MO(1) .gt. 0) then + if (trim(wfObject%name) .ne. trim(CONTROL_instance%IONIZE_SPECIES(1)) ) return + activeOrbitals = 0 + do i= 1, size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .lt. 1.0_8 ) activeOrbitals=activeOrbitals+1 + end do + call Vector_constructorInteger(orbitalsVector,activeOrbitals) + ii=0 + do i= 1, size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .lt. 1.0_8 ) then + ii=ii+1 + orbitalsVector%values(ii)=CONTROL_instance%IONIZE_MO(i) + end if + end do + end if + call Matrix_copyConstructor(auxOverlapMatrix,wfObject%overlapMatrix) - occupatedOrbitals = MolecularSystem_getOcupationNumber(wfObject%species ) - total3 = int(occupatedOrbitals,8) + numberOfContractions = MolecularSystem_getTotalnumberOfContractions(wfObject%species ) - call Matrix_constructor (matchingMatrix, total3, total3) + call Matrix_constructor (matchingMatrix, int(activeOrbitals,8), int(activeOrbitals,8)) matchingMatrix%values = & matmul( matmul( transpose( previousWavefunctionCoefficients%values ) , & @@ -421,26 +398,28 @@ subroutine SingleSCF_orbitalExchange(wfObject,previousWavefunctionCoefficients) astrayOrbitals=0 !!! number of orbitals that must be reordered - do ii= 1, occupatedOrbitals - + do i= 1, activeOrbitals + ii = orbitalsVector%values(i) + if ( abs(matchingMatrix%values(ii,ii)) < threshold ) astrayOrbitals=astrayOrbitals+1 !DEBUG - ! print *,"Antes ","Orbital: ",ii," ",abs(matchingMatrix%values(ii,ii))," energy ",wfObject%molecularOrbitalsEnergy%values(ii) - - if ( abs(matchingMatrix%values(ii,ii)) < threshold ) then - astrayOrbitals=astrayOrbitals+1 - end if - + if ( CONTROL_instance%DEBUG_SCFS) print *,"Antes ","Orbital: ",ii," ",abs(matchingMatrix%values(ii,ii))," energy ",wfObject%molecularOrbitalsEnergy%values(ii) end do !DEBUG - ! print *,"Number of astrayOrbitals",astrayOrbitals - ! print *,"Initial MM:" - ! call Matrix_show (matchingMatrix) + if (CONTROL_instance%DEBUG_SCFS) print *,"Number of astrayOrbitals",astrayOrbitals - call Vector_constructor (deltaVector,astrayOrbitals) + if(astrayOrbitals .eq. 0 ) return + if (CONTROL_instance%DEBUG_SCFS) then + print *,"Initial MM:" + call Matrix_show (matchingMatrix) + end if + + call Vector_constructor (deltaVector,astrayOrbitals) + jj=0 - do ii= 1, occupatedOrbitals + do i= 1, activeOrbitals + ii = orbitalsVector%values(i) if ( abs(matchingMatrix%values(ii,ii)) < threshold ) then jj=jj+1 deltaVector%values(jj)=ii @@ -449,61 +428,60 @@ subroutine SingleSCF_orbitalExchange(wfObject,previousWavefunctionCoefficients) prodigals = astrayOrbitals - if(astrayOrbitals .ge. 1 .and. wfObject%numberOfIterations > startIteration) then - ! print *, "Switching orbitals..." - - do ii=1,astrayOrbitals - - call Matrix_copyConstructor( auxiliaryMatrix, wfObject%waveFunctionCoefficients ) - call Matrix_copyConstructor( matchingMatrix2, matchingMatrix ) - - search=deltaVector%values(ii) - ! print *,"search: ",search - - do jj=1, numberOfContractions - - !trial=deltaVector%values(jj) - trial=jj - ! print *,"trial: ",trial + if(wfObject%numberOfIterations .lt. startIteration) return + + if (CONTROL_instance%DEBUG_SCFS) print *, "Starting switch orbitals search..." - ! print *,"valor: ",abs(matchingMatrix%values(search,trial)) + do ii=1,astrayOrbitals - if ( abs(matchingMatrix%values(search,trial)) > threshold ) then + call Matrix_copyConstructor( auxiliaryMatrix, wfObject%waveFunctionCoefficients ) + call Matrix_copyConstructor( matchingMatrix2, matchingMatrix ) - wfObject%waveFunctionCoefficients%values(:,search)=auxiliaryMatrix%values(:,trial) - wfObject%waveFunctionCoefficients%values(:,trial)=auxiliaryMatrix%values(:,search) + search=deltaVector%values(ii) + ! print *,"search: ",search + maxOverlap=0.0 + do jj=1, numberOfContractions + if( abs(matchingMatrix%values(search,jj)) > maxOverlap ) then + maxOverlap = abs(matchingMatrix%values(search,jj)) + trial = jj + end if + end do - hold=wfObject%molecularOrbitalsEnergy%values(search) - wfObject%molecularOrbitalsEnergy%values(search)=wfObject%molecularOrbitalsEnergy%values(trial) - wfObject%molecularOrbitalsEnergy%values(trial)=hold + if ( CONTROL_instance%DEBUG_SCFS) then + print *, "search", search, "trial: ",trial, "valor: ",abs(matchingMatrix%values(search,trial)) + end if - matchingMatrix%values(:,trial)=matchingMatrix2%values(:,search) - matchingMatrix%values(:,search)=matchingMatrix2%values(:,trial) + wfObject%waveFunctionCoefficients%values(:,search)=auxiliaryMatrix%values(:,trial) + wfObject%waveFunctionCoefficients%values(:,trial)=auxiliaryMatrix%values(:,search) - prodigals = prodigals - 1 + hold=wfObject%molecularOrbitalsEnergy%values(search) + wfObject%molecularOrbitalsEnergy%values(search)=wfObject%molecularOrbitalsEnergy%values(trial) + wfObject%molecularOrbitalsEnergy%values(trial)=hold - print *, "Switching orbital... ",search," with ",trial, " for ", trim(wfObject%name) + matchingMatrix%values(:,trial)=matchingMatrix2%values(:,search) + matchingMatrix%values(:,search)=matchingMatrix2%values(:,trial) - exit + prodigals = prodigals - 1 - end if + if (CONTROL_instance%PRINT_LEVEL>0) print *, "Switching orbital... ",search," with ",trial, " for ", trim(wfObject%name), " overlap", maxOverlap - end do + ! if (prodigals<2) then + ! exit + ! end if - ! if (prodigals<2) then - ! exit - ! end if + end do - end do + if (CONTROL_instance%DEBUG_SCFS) then + print *, "Coefficients after switching for ", wfObject%species + call Matrix_show(wfObject%waveFunctionCoefficients) + end if - if ( CONTROL_instance%DEBUG_SCFS) then - print *, "Coefficients after switching for ", wfObject%species - call Matrix_show(wfObject%waveFunctionCoefficients) - end if + call Matrix_destructor(matchingMatrix) + call Matrix_destructor(matchingMatrix2) + call Matrix_destructor(auxiliaryMatrix) + call Vector_destructor(deltaVector) + call Vector_destructorInteger(orbitalsVector) - call Matrix_destructor(matchingMatrix) - call Vector_destructor(deltaVector) - end if end subroutine SingleSCF_OrbitalExchange !!********************************************************************************************** @@ -511,21 +489,9 @@ end subroutine SingleSCF_OrbitalExchange subroutine SingleSCF_readFrozenOrbitals(wfObject) type(WaveFunction) :: wfObject - type(Matrix) :: fockMatrixTransformed type(Matrix) :: auxiliaryMatrix - type(Matrix) :: previousWavefunctionCoefficients, matchingMatrix, matchingMatrix2 - type(Vector) :: deltaVector - type(Matrix) :: auxOverlapMatrix - - integer(8) :: total3 integer(8) :: numberOfContractions - real(8) :: threshold, hold - real(8) :: levelShiftingFactor, normCheck - integer :: ocupados, occupatedOrbitals, prodigals - integer :: totales - integer :: search, trial - integer :: ii, jj, astrayOrbitals, startIteration - integer :: i,j, mu, nu, index + integer :: i,mu, nu logical :: existFile character(50) :: wfnFile diff --git a/src/scf/WaveFunction.f90 b/src/scf/WaveFunction.f90 index d345874d..6f04e641 100644 --- a/src/scf/WaveFunction.f90 +++ b/src/scf/WaveFunction.f90 @@ -1760,6 +1760,7 @@ subroutine WaveFunction_buildDensityMatrix(this) implicit none type(WaveFunction) :: this + type(Matrix) :: auxMatrix integer :: orderMatrix integer :: ocupationNumber integer :: i @@ -1772,33 +1773,31 @@ subroutine WaveFunction_buildDensityMatrix(this) this%densityMatrix%values = 0.0_8 - !! Segment for fractional occupations: 1 - if (CONTROL_instance%IONIZE_MO /= 0 .and. trim(this%name) == trim(CONTROL_instance%IONIZE_SPECIE(1)) ) then - this%waveFunctionCoefficients%values(:,CONTROL_instance%IONIZE_MO) = & - this%waveFunctionCoefficients%values(:,CONTROL_instance%IONIZE_MO)*sqrt(CONTROL_instance%MO_FRACTION_OCCUPATION) - + call Matrix_copyConstructor(auxMatrix, this%waveFunctionCoefficients) + + !! Segment for fractional occupations: introduce fractional occupation + if (trim(this%name) == trim(CONTROL_instance%IONIZE_SPECIES(1)) ) then + do i=1,size(CONTROL_instance%IONIZE_MO) + if(CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .lt. 1.0_8) & + auxMatrix%values(:,CONTROL_instance%IONIZE_MO(i)) = & + auxMatrix%values(:,CONTROL_instance%IONIZE_MO(i))*sqrt(CONTROL_instance%MO_FRACTION_OCCUPATION(i)) + end do end if - + do i = 1 , orderMatrix do j = 1 , orderMatrix do k = 1 , ocupationNumber this%densityMatrix%values(i,j) = & this%densityMatrix%values( i,j ) + & - ( this%waveFunctionCoefficients%values(i,k) & - * this%waveFunctionCoefficients%values(j,k) ) + ( auxMatrix%values(i,k) & + * auxMatrix%values(j,k) ) end do end do end do this%densityMatrix%values = MolecularSystem_getEta(this%species) * this%densityMatrix%values - !! Segment for fractional occupations: 1 - if (CONTROL_instance%IONIZE_MO /= 0 .and. trim(this%name) == trim(CONTROL_instance%IONIZE_SPECIE(1))) then - this%waveFunctionCoefficients%values(:,CONTROL_instance%IONIZE_MO) = & - this%waveFunctionCoefficients%values(:,CONTROL_instance%IONIZE_MO)/sqrt(CONTROL_instance%MO_FRACTION_OCCUPATION) - end if - !!DEBUG if ( CONTROL_instance%DEBUG_SCFS) then print *,"Density Matrix ", trim(this%name) @@ -2628,6 +2627,7 @@ subroutine WaveFunction_cosmoQuantumCharge() write(*,*) "COSMO Charges for ",MolecularSystem_getNameOfSpecies( f )," = ", sum(qiCosmo(:)) end do + close(wfnUnit) charges_file="qTotalCosmo.charges" @@ -2639,7 +2639,47 @@ subroutine WaveFunction_cosmoQuantumCharge() ! write(*,*)"sumqTotalCosmo ", sum(qTotalCosmo(:)) end subroutine WaveFunction_cosmoQuantumCharge + + !> + !! @brief Remove orbitals from basis functions deleted by the overlap criteria selected + subroutine Wavefunction_removeOrbitalsBelowEigenThreshold(this) + implicit none + type(WaveFunction) :: this + + integer(8) :: numberOfContractions + real(8) :: normCheck + integer :: i, j, mu, nu, index + + numberOfContractions = MolecularSystem_getTotalnumberOfContractions(this%species) + + i=0 + do index = 1 , numberOfContractions + i=i+1 + normCheck=0.0 + if ( abs(this%molecularOrbitalsEnergy%values(i)) .lt. CONTROL_instance%OVERLAP_EIGEN_THRESHOLD ) then + do mu = 1 , numberOfContractions + do nu = 1 , numberOfContractions + normCheck=normCheck+this%waveFunctionCoefficients%values(mu,i)*& + this%waveFunctionCoefficients%values(nu,i)*& + this%overlapMatrix%values(mu,nu) + end do + end do + ! print *, "eigenvalue", i, this%molecularOrbitalsEnergy%values(i), "normCheck", normCheck + if ( normCheck .lt. CONTROL_instance%OVERLAP_EIGEN_THRESHOLD) then + ! Shift orbital coefficients to the end of the matrix and Make energy a very large number + do j = i , numberOfContractions-1 + this%molecularOrbitalsEnergy%values(j)=this%molecularOrbitalsEnergy%values(j+1) + this%waveFunctionCoefficients%values(:,j) = this%waveFunctionCoefficients%values(:,j+1) + end do + this%molecularOrbitalsEnergy%values(numberOfContractions)=1/CONTROL_instance%OVERLAP_EIGEN_THRESHOLD + this%waveFunctionCoefficients%values(:,numberOfContractions)=0.0 + i=i-1 + end if + + end if + end do + end subroutine Wavefunction_removeOrbitalsBelowEigenThreshold end module WaveFunction_ diff --git a/test/H2O.APMO.UP2.IS.IMO1.IT-E.lowdin b/test/H2O.APMO.UP2.IS.IMO1.IT-E.lowdin index f07a6094..f5216d98 100644 --- a/test/H2O.APMO.UP2.IS.IMO1.IT-E.lowdin +++ b/test/H2O.APMO.UP2.IS.IMO1.IT-E.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="H-A_1" "H-B_1" +IonizeSpecies="H-A_1" "H-B_1" IonizeMO = 1 integralsTransformationMethod = "E" END CONTROL diff --git a/test/H2O.APMO.UP2.IS.IMO1.lowdin b/test/H2O.APMO.UP2.IS.IMO1.lowdin index c4a78eb8..65481da6 100644 --- a/test/H2O.APMO.UP2.IS.IMO1.lowdin +++ b/test/H2O.APMO.UP2.IS.IMO1.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="H-A_1" "H-B_1" +IonizeSpecies="H-A_1" "H-B_1" IonizeMO = 1 END CONTROL diff --git a/test/H2O.APMO.UP2.IS.IMO2.IT-E.lowdin b/test/H2O.APMO.UP2.IS.IMO2.IT-E.lowdin index 775522ca..e61c5af3 100644 --- a/test/H2O.APMO.UP2.IS.IMO2.IT-E.lowdin +++ b/test/H2O.APMO.UP2.IS.IMO2.IT-E.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="H-A_1" "H-B_1" +IonizeSpecies="H-A_1" "H-B_1" IonizeMO = 2 integralsTransformationMethod = "E" END CONTROL diff --git a/test/H2O.APMO.UP2.IS.IMO2.lowdin b/test/H2O.APMO.UP2.IS.IMO2.lowdin index df983759..e6414160 100644 --- a/test/H2O.APMO.UP2.IS.IMO2.lowdin +++ b/test/H2O.APMO.UP2.IS.IMO2.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="H-A_1" "H-B_1" +IonizeSpecies="H-A_1" "H-B_1" IonizeMO = 2 END CONTROL diff --git a/test/H2O.APMO.UP2.IS.lowdin b/test/H2O.APMO.UP2.IS.lowdin index e784f6d2..ab25012f 100644 --- a/test/H2O.APMO.UP2.IS.lowdin +++ b/test/H2O.APMO.UP2.IS.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="H-A_1" "H-B_1" +IonizeSpecies="H-A_1" "H-B_1" END CONTROL diff --git a/test/H2O.APMO.UP2.IS2.IT-E.lowdin b/test/H2O.APMO.UP2.IS2.IT-E.lowdin index 41857455..2c66bff6 100644 --- a/test/H2O.APMO.UP2.IS2.IT-E.lowdin +++ b/test/H2O.APMO.UP2.IS2.IT-E.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="E-ALPHA" "E-BETA" "H-A_1" +IonizeSpecies="E-ALPHA" "E-BETA" "H-A_1" integralsTransformationMethod = "E" END CONTROL diff --git a/test/H2O.APMO.UP2.IS2.lowdin b/test/H2O.APMO.UP2.IS2.lowdin index 9b60f51c..2555850a 100644 --- a/test/H2O.APMO.UP2.IS2.lowdin +++ b/test/H2O.APMO.UP2.IS2.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -IonizeSpecie="E-ALPHA" "E-BETA" "H-A_1" +IonizeSpecies="E-ALPHA" "E-BETA" "H-A_1" END CONTROL diff --git a/test/H2O.BOA.UP2.IT-E.lowdin b/test/H2O.BOA.UP2.IT-E.lowdin index a6b1b501..4f320df7 100644 --- a/test/H2O.BOA.UP2.IT-E.lowdin +++ b/test/H2O.BOA.UP2.IT-E.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -! ionizeSpecie="E-ALPHA" +! ionizeSpecies="E-ALPHA" ! ionizeMO=2 ! ptJustOneOrbital=T ! ptTransitionOperator=T diff --git a/test/H2O.BOA.UP2.lowdin b/test/H2O.BOA.UP2.lowdin index 11a81466..a513b118 100644 --- a/test/H2O.BOA.UP2.lowdin +++ b/test/H2O.BOA.UP2.lowdin @@ -29,7 +29,7 @@ END TASKS CONTROL readCoefficients=F -! ionizeSpecie="E-ALPHA" +! ionizeSpecies="E-ALPHA" ! ionizeMO=2 ! ptJustOneOrbital=T ! ptTransitionOperator=T diff --git a/test/H2O.coreHole.lowdin b/test/H2O.coreHole.lowdin new file mode 100644 index 00000000..0e84c262 --- /dev/null +++ b/test/H2O.coreHole.lowdin @@ -0,0 +1,29 @@ +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Test RHF +% H2O Single Point +% Electronic Hartree Fock +% Basis sets used: +% 6-311g e- +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SYSTEM_DESCRIPTION='Molecula de H2O' + +GEOMETRY + e-[O] 6-311G 0.000000 0.000000 -0.066575 + e-[H] 6-311G 0.000000 0.754175 0.528381 + e-[H] 6-311G 0.000000 -0.754174 0.528382 + O dirac 0.000000 0.000000 -0.066575 + H_1 dirac 0.000000 0.754175 0.528381 + H_1 dirac 0.000000 -0.754174 0.528382 +END GEOMETRY + +TASKS + method = "UHF" +END TASKS + +CONTROL +readCoefficients=F +ionizeMO=1 +ionizeSpecies="E-BETA" +MOfractionOccupation=0.0 +END CONTROL diff --git a/test/H2O.coreHole.py b/test/H2O.coreHole.py new file mode 100644 index 00000000..30f34846 --- /dev/null +++ b/test/H2O.coreHole.py @@ -0,0 +1,46 @@ +#!/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 = "H2O.coreHole" +inputName = testName + ".lowdin" +outputName = testName + ".out" + +# Reference values + +refTotalEnergy = -56.177665367588 + +# 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 + +for line in outputRead: + if "TOTAL ENERGY =" in line: + totalEnergy = float(line.split()[3]) + +diffTotalEnergy = abs(refTotalEnergy - totalEnergy) + +if (diffTotalEnergy <= 1E-8): + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + print("Difference HF: " + str(diffTotalEnergy)) + sys.exit(1) + +output.close() diff --git a/test/HCN.APMO.RKS.parallel.lowdin b/test/HCN.APMO.RKS.parallel.lowdin new file mode 100644 index 00000000..f02fa3be --- /dev/null +++ b/test/HCN.APMO.RKS.parallel.lowdin @@ -0,0 +1,35 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Test RKS +% HCN Single Point +% APMO DFT with nuclear electron correlation +% Results +% refExchangeCorrelationEnergy=-10.15758267 +% refNuclearElectronCorrelationEnergy=-0.03005334 +% refTotalEnergy=-93.367393799928 +% Author: Felix Moncada - QCC-UNAL/2018 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SYSTEM_DESCRIPTION='Molecula de H2' + +GEOMETRY +e-(N) CC-PVDZ 0.0000 0.00000 1.156 +e-(C) CC-PVDZ 0.0000 0.00000 0.0000 +e-(H) CC-PVDZ 0.0000 0.00000 -1.064 +N dirac 0.0000 0.00000 1.156 +C dirac 0.0000 0.00000 0.0000 +H_1 Nakai-3-SP 0.0000 0.00000 -1.064 +END GEOMETRY + +TASKS +method = "RKS" +END TASKS + +CONTROL +readCoefficients=.F. +electronExchangeCorrelationFunctional="B3LYP" +nuclearElectronCorrelationFunctional="epc17-2" +gridStorage="MEMORY" +END CONTROL + + + diff --git a/test/HCOOPs.HF.DensCube.lowdin b/test/HCOOPs.HF.DensCube.lowdin index fd3267ac..0626cb09 100644 --- a/test/HCOOPs.HF.DensCube.lowdin +++ b/test/HCOOPs.HF.DensCube.lowdin @@ -21,7 +21,8 @@ CONTROL END CONTROL OUTPUTS - densityCube cubeSize=20 point1=0.0 0.0 -2.5 + densityCube cubeSize=5 point1=0.0 0.0 0.0 species="E-" + densityCube cubeSize=20 point1=0.0 0.0 -2.5 species="POSITRON" END OUTPUTS diff --git a/test/HCOOPs.HF.DensCube.py b/test/HCOOPs.HF.DensCube.py index f1a6b6d1..43ac5465 100644 --- a/test/HCOOPs.HF.DensCube.py +++ b/test/HCOOPs.HF.DensCube.py @@ -13,13 +13,13 @@ inputName = testName + ".lowdin" outputName = testName + ".out" cube1Name = testName + ".E-.dens.cub" -cube2Name = testName + ".POSITRON.dens.cub" +cube2Name = testName + ".POSITRON.2.dens.cub" # Reference values and tolerance refValues = { "HF energy" : [-188.362545831570,1E-8], -"Num e- in cube" : [23.98744049,1E-1], -"Num e+ in cube" : [0.96711581,1E-2], +"Num e- in cube" : [24.0,1E-1], +"Num e+ in cube" : [1.0,1E-2], } testValues = dict(refValues) #copy diff --git a/test/LiF2-.EDA.lowdin b/test/LiF2-.EDA.lowdin new file mode 100644 index 00000000..f58da1f1 --- /dev/null +++ b/test/LiF2-.EDA.lowdin @@ -0,0 +1,21 @@ + +GEOMETRY +e-(F) aug-cc-pVDZ 0.0000 0.0000 -1.70 addParticles=1 +e-(Li) aug-cc-pVDZ 0.0000 0.0000 0.0 +e-(F) aug-cc-pVDZ 0.0000 0.0000 1.70 +e+ PSX-TZ 0.0000 0.0000 -1.70 addParticles=-2 +e+ PSX-TZ 0.0000 0.0000 1.70 +F dirac 0.0000 0.0000 -1.70 +Li dirac 0.0000 0.0000 0.0 +F dirac 0.0000 0.0000 1.70 +END GEOMETRY + +TASKS + method = "RHF" + mollerPlessetCorrection=2 +END TASKS + +CONTROL + readCoefficients=.F. + scfTotalEnergyTolerance=1E-12 +END CONTROL diff --git a/test/LiF2-.MP2.lowdin b/test/LiF2-.MP2.lowdin new file mode 100644 index 00000000..34921b23 --- /dev/null +++ b/test/LiF2-.MP2.lowdin @@ -0,0 +1,21 @@ + +GEOMETRY +e-(F) aug-cc-pVDZ 0.0000 0.0000 -1.70 addParticles=1 +e-(Li) aug-cc-pVDZ 0.0000 0.0000 0.0 +e-(F) aug-cc-pVDZ 0.0000 0.0000 1.70 +e+ PSX-TZ 0.0000 0.0000 -1.70 addParticles=-2 +e+ PSX-TZ 0.0000 0.0000 1.70 +F dirac 0.0000 0.0000 -1.70 +Li dirac 0.0000 0.0000 0.0 +F dirac 0.0000 0.0000 1.70 +END GEOMETRY + +TASKS + method = "RHF" + mollerPlessetCorrection=2 +END TASKS + +CONTROL + readCoefficients=.F. + totalEnergyTolerance=1E-12 +END CONTROL diff --git a/test/LiF2-.addParticles.lowdin b/test/LiF2-.addParticles.lowdin new file mode 100644 index 00000000..0a9f13c2 --- /dev/null +++ b/test/LiF2-.addParticles.lowdin @@ -0,0 +1,20 @@ + +GEOMETRY +e-(F) aug-cc-pVDZ 0.0000 0.0000 -1.70 addParticles=1 +e-(Li) aug-cc-pVDZ 0.0000 0.0000 0.0 +e-(F) aug-cc-pVDZ 0.0000 0.0000 1.70 +e+ PSX-TZ 0.0000 0.0000 -1.70 addParticles=-1 +e+ PSX-TZ 0.0000 0.0000 1.70 addParticles=-1 +F dirac 0.0000 0.0000 -1.70 +Li dirac 0.0000 0.0000 0.0 +F dirac 0.0000 0.0000 1.70 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=.F. + totalEnergyTolerance=1E-12 +END CONTROL diff --git a/test/LiF2-e+.EDA.lowdin b/test/LiF2-e+.EDA.lowdin new file mode 100644 index 00000000..7a1069d2 --- /dev/null +++ b/test/LiF2-e+.EDA.lowdin @@ -0,0 +1,21 @@ + +GEOMETRY +e-(F) aug-cc-pVDZ 0.0000 0.0000 -1.70 addParticles=1 +e-(Li) aug-cc-pVDZ 0.0000 0.0000 0.0 +e-(F) aug-cc-pVDZ 0.0000 0.0000 1.70 +e+ PSX-TZ 0.0000 0.0000 -1.70 addParticles=-1 +e+ PSX-TZ 0.0000 0.0000 1.70 +F dirac 0.0000 0.0000 -1.70 +Li dirac 0.0000 0.0000 0.0 +F dirac 0.0000 0.0000 1.70 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=.T. + noSCF=.T. + scfTotalEnergyTolerance=1E-12 +END CONTROL diff --git a/test/LiF2-e+.MP2.lowdin b/test/LiF2-e+.MP2.lowdin new file mode 100644 index 00000000..41c2f65e --- /dev/null +++ b/test/LiF2-e+.MP2.lowdin @@ -0,0 +1,21 @@ + +GEOMETRY +e-(F) aug-cc-pVDZ 0.0000 0.0000 -1.70 addParticles=1 +e-(Li) aug-cc-pVDZ 0.0000 0.0000 0.0 +e-(F) aug-cc-pVDZ 0.0000 0.0000 1.70 +e+ PSX-TZ 0.0000 0.0000 -1.70 addParticles=-1 +e+ PSX-TZ 0.0000 0.0000 1.70 +F dirac 0.0000 0.0000 -1.70 +Li dirac 0.0000 0.0000 0.0 +F dirac 0.0000 0.0000 1.70 +END GEOMETRY + +TASKS + method = "RHF" + mollerPlessetCorrection=2 +END TASKS + +CONTROL + readCoefficients=.F. + totalEnergyTolerance=1E-12 +END CONTROL diff --git a/test/LiF2-e+.ghost.lowdin b/test/LiF2-e+.ghost.lowdin new file mode 100644 index 00000000..1e5a1ae6 --- /dev/null +++ b/test/LiF2-e+.ghost.lowdin @@ -0,0 +1,21 @@ + +GEOMETRY +e-(F) aug-cc-pVDZ 0.0000 0.0000 -1.70 addParticles=1 +e-(Li) aug-cc-pVDZ 0.0000 0.0000 0.0 +e-(F) aug-cc-pVDZ 0.0000 0.0000 1.70 +e+ PSX-TZ 0.0000 0.0000 -1.70 addParticles=-1 +e+ PSX-TZ 0.0000 0.0000 1.70 +F dirac 0.0000 0.0000 -1.70 +Li dirac 0.0000 0.0000 0.0 +F dirac 0.0000 0.0000 1.70 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=.F. + scfGhostSpecies="POSITRON" + TotalEnergyTolerance=1E-12 +END CONTROL diff --git a/test/LiH.APMO.TOP2.IT-E.lowdin b/test/LiH.APMO.TOP2.IT-E.lowdin index 64cca285..4e5311da 100644 --- a/test/LiH.APMO.TOP2.IT-E.lowdin +++ b/test/LiH.APMO.TOP2.IT-E.lowdin @@ -16,7 +16,7 @@ CONTROL readCoefficients=F ! active to compare with Romero.JCP.137.74105.2012 ! buildTwoParticlesMatrixForOneParticle=T - ionizeSpecie="E-ALPHA" + ionizeSpecies="E-ALPHA" IonizeMO=2 MOfractionOccupation=0.5 ptTransitionOperator=T diff --git a/test/LiH.APMO.TOP2.lowdin b/test/LiH.APMO.TOP2.lowdin index 03153f5a..f22eacfd 100644 --- a/test/LiH.APMO.TOP2.lowdin +++ b/test/LiH.APMO.TOP2.lowdin @@ -16,7 +16,7 @@ CONTROL readCoefficients=F ! active to compare with Romero.JCP.137.74105.2012 ! buildTwoParticlesMatrixForOneParticle=T - ionizeSpecie="E-ALPHA" + ionizeSpecies="E-ALPHA" IonizeMO=2 MOfractionOccupation=0.5 ptTransitionOperator=T diff --git a/test/LiH.TOP2.IT-E.lowdin b/test/LiH.TOP2.IT-E.lowdin index 80f5078a..82bc9c29 100644 --- a/test/LiH.TOP2.IT-E.lowdin +++ b/test/LiH.TOP2.IT-E.lowdin @@ -15,7 +15,7 @@ END TASKS CONTROL readCoefficients=F !buildTwoParticlesMatrixForOneParticle=T ! active to compare with Romero.JCP.137.74105.2012 - ionizeSpecie="E-ALPHA" + ionizeSpecies="E-ALPHA" IonizeMO=2 MOfractionOccupation=0.5 ptTransitionOperator=T diff --git a/test/LiH.TOP2.lowdin b/test/LiH.TOP2.lowdin index cefa412d..38b68065 100644 --- a/test/LiH.TOP2.lowdin +++ b/test/LiH.TOP2.lowdin @@ -15,7 +15,7 @@ END TASKS CONTROL readCoefficients=F !buildTwoParticlesMatrixForOneParticle=T ! active to compare with Romero.JCP.137.74105.2012 - ionizeSpecie="E-ALPHA" + ionizeSpecies="E-ALPHA" IonizeMO=2 MOfractionOccupation=0.5 ptTransitionOperator=T diff --git a/test/Mg.e+.gribakin.lowdin b/test/Mg.e+.gribakin.lowdin new file mode 100644 index 00000000..e3eba0d6 --- /dev/null +++ b/test/Mg.e+.gribakin.lowdin @@ -0,0 +1,26 @@ + +GEOMETRY + e-(Mg) 6-311++G.d.p 0 0 0 + Mg dirac 0 0 0 + E+ GRIBAKIN-10S4P2D 0 0 0.00 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +EXTERPOTENTIAL + E- GRIBAKIN0 + E+ GRIBAKINMG +END EXTERPOTENTIAL + +CONTROL + units="BOHRS" +readCoefficients=F +convergenceMethod=1 + overlapEigenThreshold=1E-6 + SCFGhostSpecies="POSITRON" + HFprintEigenvalues=T + totalEnergyTolerance = 1E-9 +END CONTROL + diff --git a/test/Mg2.e+.gribakin.lowdin b/test/Mg2.e+.gribakin.lowdin new file mode 100644 index 00000000..d5ea4c1e --- /dev/null +++ b/test/Mg2.e+.gribakin.lowdin @@ -0,0 +1,29 @@ + +GEOMETRY + e-(Mg) 6-311++G.d.p 0 0 0 + Mg dirac 0 0 0 + e-(Mg) 6-311++G.d.p 0 0 6.001 + Mg dirac 0 0 6.001 + E+ GRIBAKIN-10S4P2D 0 0 0.00 + E+ GRIBAKIN-10S4P2D 0 0 6.001 addParticles=-1 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +EXTERPOTENTIAL + E- GRIBAKIN0 + E+ GRIBAKINMG2 +END EXTERPOTENTIAL + +CONTROL + units="BOHRS" +readCoefficients=F +convergenceMethod=1 + overlapEigenThreshold=1E-6 + SCFGhostSpecies="POSITRON" + HFprintEigenvalues=T + totalEnergyTolerance = 1E-9 +END CONTROL + diff --git a/test/PsF-RENP3.lowdin b/test/PsF-RENP3.lowdin index 2e9cb3ab..c47cbe70 100644 --- a/test/PsF-RENP3.lowdin +++ b/test/PsF-RENP3.lowdin @@ -11,7 +11,7 @@ END TASKS CONTROL ptP3Method = "Ren-P3" - ionizeSpecie = "positron" + ionizeSpecies = "positron" ionizeMO= 1 ptJustOneOrbital=.T. END CONTROL diff --git a/test/PsH.FCI.test.lowdin b/test/PsH.FCI.test.lowdin new file mode 100644 index 00000000..016bd225 --- /dev/null +++ b/test/PsH.FCI.test.lowdin @@ -0,0 +1,33 @@ +GEOMETRY + e-(H) SHARON-E-6S2P 0.00 0.00 0.00 addParticles=1 + H dirac 0.00 0.00 0.00 + e+ SHARON-E+6S2P 0.00 0.00 0.00 +END GEOMETRY + +TASKS + method = "UHF" + configurationInteractionLevel ="FCI" + !configurationInteractionLevel ="CISD" +END TASKS + +CONTROL +readCoefficients=F +numberOfCIstates=3 +CINaturalOrbitals=T + CIStatesToPrint = 1 + !CIdiagonalizationMethod = "DSYEVX" + CIdiagonalizationMethod = "JADAMILU" + !CIdiagonalizationMethod = "ARPACK" + !CIPrintEigenVectorsFormat = "NONE" + !CIPrintEigenVectorsFormat = "OCCUPIED" + CIPrintEigenVectorsFormat = "ORBITALS" + CIPrintThreshold = 5e-2 + buildTwoParticlesMatrixForOneParticle=T +END CONTROL + +INPUT_CI + species="E-ALPHA" core=0 active=0 + species="E-BETA" core=0 active=0 + species="POSITRON" core=0 active=0 +END INPUT_CI +