diff --git a/README.md b/README.md index ba7e37f0..140b7ee0 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,6 @@ Installation notes. * A standard FORTRAN compiler. gfortran and intel FORTRAN compiler have been tested. * Lapack or MKL libraries. -* Arpack library. * LIBINT library version 2.0, and version 1.1.5 * LIBXC library @@ -63,7 +62,7 @@ To visualize the html documentation use: ### Step-by-step installation example: (replace apt-get with your preferred package manager) ### sudo apt-get update - sudo apt-get -y install wget git build-essential liblapack-dev libblas-dev libgsl0-dev autotools-dev automake libtool gfortran python3 gawk libeigen3-dev libgmp-dev libboost-all-dev libarpack2-dev + sudo apt-get -y install wget git build-essential liblapack-dev libblas-dev libgsl0-dev autotools-dev automake libtool gfortran python3 gawk libeigen3-dev libgmp-dev libboost-all-dev # Define ENV Variables export WORKDIR=$PWD/dependencies export PATH=$PATH:$WORKDIR/bin @@ -106,7 +105,7 @@ To visualize the html documentation use: cd .. # Configure Lowdin - ./configure -p $WORKDIR/bin -s /tmp -l "-lblas -llapack -larpack" + ./configure -p $WORKDIR/bin -s /tmp -l "-lblas -llapack" # Build Lowdin make -j 4 # Install Lowdin diff --git a/configure b/configure index 86fb684c..26a1a9d7 100755 --- a/configure +++ b/configure @@ -46,7 +46,7 @@ CFLAGS="-g -O2 -fopenmp -D gfortran" CPPFLAGS="-g -O2 -fopenmp -D gfortran -std=c++11 " FCFLAGS="-O2 -ftree-vectorize -ffast-math -fstack-arrays -fpic -fopenmp -cpp -ffree-line-length-none -fno-range-check -fallow-invalid-boz -fbacktrace -g -D gfortran" GSL_LIBS="-lgsl -lgslcblas" -LAPACK_LIBS="-lblas -llapack -larpack" +LAPACK_LIBS="-lblas -llapack" STDCPP_LIBS="-lstdc++" LIBINT_LIBS="-lr12 -lderiv -lint -lint2" LIBXC_LIBS="-lxcf03 -lxc" @@ -211,9 +211,9 @@ else ############################################################################### if [ "$FC" = "ifort" ] ; then - LAPACK_LIBS="-mkl -larpack" + LAPACK_LIBS="-mkl" else - LAPACK_LIBS="-lblas -llapack -larpack" + LAPACK_LIBS="-lblas -llapack" fi if [ "$CUDAFLAGS" = "yes" ] ; then diff --git a/lib/basis/AUG-CC-PVQZ b/lib/basis/AUG-CC-PVQZ index 3eabdfe9..e404dd35 100644 --- a/lib/basis/AUG-CC-PVQZ +++ b/lib/basis/AUG-CC-PVQZ @@ -1750,3 +1750,193 @@ O-SCANDIUM SC (AUG-CC-PVQZ) BASIS TYPE: 1 0.25730000 1.00000000 30 4 1 0.07236000 1.00000000 + +O-BROMINE BR (aug-cc-pVQZ) BASIS TYPE: 1 +# (88s,52p,15d,3f,2g) -> [8s,7p,5d,3f,2g] + 25 + 1 0 21 + 1.647500E+07 3.400000E-06 + 2.466600E+06 2.670000E-05 + 5.613100E+05 1.404000E-04 + 1.589900E+05 5.927000E-04 + 5.186900E+04 2.156100E-03 + 1.872600E+04 6.995900E-03 + 7.303600E+03 2.056450E-02 + 3.029100E+03 5.458930E-02 + 1.320800E+03 1.275226E-01 + 6.000300E+02 2.459780E-01 + 2.819000E+02 3.436508E-01 + 1.355400E+02 2.702530E-01 + 6.487000E+01 7.447950E-02 + 3.212900E+01 8.787000E-04 + 1.603700E+01 1.575500E-03 + 7.784900E+00 -7.602000E-04 + 3.724700E+00 3.211000E-04 + 1.758300E+00 -1.586000E-04 + 5.833100E-01 6.590000E-05 + 2.785600E-01 -3.740000E-05 + 1.182900E-01 1.000000E-05 + 2 0 21 + 1.647500E+07 -1.100000E-06 + 2.466600E+06 -8.400000E-06 + 5.613100E+05 -4.410000E-05 + 1.589900E+05 -1.862000E-04 + 5.186900E+04 -6.783000E-04 + 1.872600E+04 -2.212200E-03 + 7.303600E+03 -6.575200E-03 + 3.029100E+03 -1.793280E-02 + 1.320800E+03 -4.433210E-02 + 6.000300E+02 -9.634780E-02 + 2.819000E+02 -1.696814E-01 + 1.355400E+02 -1.920769E-01 + 6.487000E+01 2.087310E-02 + 3.212900E+01 4.744996E-01 + 1.603700E+01 5.214907E-01 + 7.784900E+00 1.348001E-01 + 3.724700E+00 3.661400E-03 + 1.758300E+00 1.884000E-03 + 5.833100E-01 -4.582000E-04 + 2.785600E-01 2.165000E-04 + 1.182900E-01 -6.780000E-05 + 3 0 21 + 1.647500E+07 4.000000E-07 + 2.466600E+06 3.300000E-06 + 5.613100E+05 1.750000E-05 + 1.589900E+05 7.400000E-05 + 5.186900E+04 2.697000E-04 + 1.872600E+04 8.799000E-04 + 7.303600E+03 2.619800E-03 + 3.029100E+03 7.167100E-03 + 1.320800E+03 1.785610E-02 + 6.000300E+02 3.939600E-02 + 2.819000E+02 7.171020E-02 + 1.355400E+02 8.588770E-02 + 6.487000E+01 -1.038610E-02 + 3.212900E+01 -3.040135E-01 + 1.603700E+01 -4.933178E-01 + 7.784900E+00 1.608900E-02 + 3.724700E+00 7.146686E-01 + 1.758300E+00 4.904795E-01 + 5.833100E-01 3.375730E-02 + 2.785600E-01 -6.792500E-03 + 1.182900E-01 1.691600E-03 + 4 0 1 + 5.833100E-01 1.000000E+00 + 5 0 1 + 2.785600E-01 1.000000E+00 + 6 0 21 + 1.647500E+07 -1.000000E-07 + 2.466600E+06 -1.000000E-06 + 5.613100E+05 -5.400000E-06 + 1.589900E+05 -2.270000E-05 + 5.186900E+04 -8.270000E-05 + 1.872600E+04 -2.694000E-04 + 7.303600E+03 -8.042000E-04 + 3.029100E+03 -2.194900E-03 + 1.320800E+03 -5.493900E-03 + 6.000300E+02 -1.209600E-02 + 2.819000E+02 -2.226230E-02 + 1.355400E+02 -2.660630E-02 + 6.487000E+01 2.758000E-03 + 3.212900E+01 1.016803E-01 + 1.603700E+01 1.704132E-01 + 7.784900E+00 -6.222000E-03 + 3.724700E+00 -3.452570E-01 + 1.758300E+00 -4.234840E-01 + 5.833100E-01 3.354720E-01 + 2.785600E-01 6.779724E-01 + 1.182900E-01 2.523080E-01 + 7 0 1 + 1.182900E-01 1.000000E+00 + 8 0 1 + 0.0442700 1.0000000 + 2 1 16 + 2.660700E+04 6.190000E-05 + 6.298200E+03 5.499000E-04 + 2.045200E+03 3.162000E-03 + 7.821600E+02 1.379790E-02 + 3.316300E+02 4.798120E-02 + 1.511100E+02 1.315710E-01 + 7.239200E+01 2.685861E-01 + 3.586200E+01 3.683473E-01 + 1.813400E+01 2.711363E-01 + 9.043000E+00 7.622220E-02 + 4.450000E+00 4.674900E-03 + 2.166100E+00 1.256500E-03 + 9.962800E-01 -2.357000E-04 + 4.544300E-01 1.332000E-04 + 1.940400E-01 -6.890000E-05 + 7.899700E-02 1.380000E-05 + 3 1 16 + 2.660700E+04 -2.480000E-05 + 6.298200E+03 -2.212000E-04 + 2.045200E+03 -1.273600E-03 + 7.821600E+02 -5.617900E-03 + 3.316300E+02 -1.986000E-02 + 1.511100E+02 -5.655310E-02 + 7.239200E+01 -1.209479E-01 + 3.586200E+01 -1.773098E-01 + 1.813400E+01 -9.214720E-02 + 9.043000E+00 2.187683E-01 + 4.450000E+00 4.854670E-01 + 2.166100E+00 3.721970E-01 + 9.962800E-01 7.769070E-02 + 4.544300E-01 6.541000E-04 + 1.940400E-01 1.799800E-03 + 7.899700E-02 -6.550000E-05 + 4 1 1 + 4.544300E-01 1.000000E+00 + 5 1 16 + 2.660700E+04 6.400000E-06 + 6.298200E+03 5.720000E-05 + 2.045200E+03 3.297000E-04 + 7.821600E+02 1.456200E-03 + 3.316300E+02 5.159100E-03 + 1.511100E+02 1.476170E-02 + 7.239200E+01 3.176940E-02 + 3.586200E+01 4.706800E-02 + 1.813400E+01 2.238710E-02 + 9.043000E+00 -7.202540E-02 + 4.450000E+00 -1.626429E-01 + 2.166100E+00 -1.496503E-01 + 9.962800E-01 1.064517E-01 + 4.544300E-01 4.360266E-01 + 1.940400E-01 4.680523E-01 + 7.899700E-02 1.691513E-01 + 6 1 1 + 1.940400E-01 1.000000E+00 + 7 1 1 + 7.899700E-02 1.000000E+00 + 8 1 1 + 0.0305130 1.0000000 + 3 2 11 + 1.289600E+03 1.190000E-04 + 3.897500E+02 1.155100E-03 + 1.517600E+02 6.764800E-03 + 6.722300E+01 2.730170E-02 + 3.191300E+01 8.092980E-02 + 1.585700E+01 1.794011E-01 + 8.054500E+00 2.840086E-01 + 4.088700E+00 3.266797E-01 + 2.055600E+00 2.584900E-01 + 9.950900E-01 1.128946E-01 + 4.231300E-01 1.726050E-02 + 4 2 1 + 9.950900E-01 1.000000E+00 + 5 2 1 + 4.231300E-01 1.000000E+00 + 6 2 1 + 1.779000E-01 1.000000E+00 + 7 2 1 + 0.0829000 1.0000000 + 4 3 1 + 8.257000E-01 1.000000E+00 + 5 3 1 + 3.407000E-01 1.000000E+00 + 6 3 1 + 0.1748000 1.0000000 + 5 4 1 + 6.491000E-01 1.0000000 + 6 4 1 + 0.3110000 1.0000000 + diff --git a/lib/basis/E+-BR-AUG-CC-PVQZ b/lib/basis/E+-BR-AUG-CC-PVQZ new file mode 100644 index 00000000..dd3c6b9f --- /dev/null +++ b/lib/basis/E+-BR-AUG-CC-PVQZ @@ -0,0 +1,189 @@ +O-POSITRON E+ (BR-AUG-CC-PVQZ) BASIS TYPE:2 +# (88s,52p,15d,3f,2g) -> [8s,7p,5d,3f,2g] + 25 + 1 0 21 + 1.647500E+07 3.400000E-06 + 2.466600E+06 2.670000E-05 + 5.613100E+05 1.404000E-04 + 1.589900E+05 5.927000E-04 + 5.186900E+04 2.156100E-03 + 1.872600E+04 6.995900E-03 + 7.303600E+03 2.056450E-02 + 3.029100E+03 5.458930E-02 + 1.320800E+03 1.275226E-01 + 6.000300E+02 2.459780E-01 + 2.819000E+02 3.436508E-01 + 1.355400E+02 2.702530E-01 + 6.487000E+01 7.447950E-02 + 3.212900E+01 8.787000E-04 + 1.603700E+01 1.575500E-03 + 7.784900E+00 -7.602000E-04 + 3.724700E+00 3.211000E-04 + 1.758300E+00 -1.586000E-04 + 5.833100E-01 6.590000E-05 + 2.785600E-01 -3.740000E-05 + 1.182900E-01 1.000000E-05 + 2 0 21 + 1.647500E+07 -1.100000E-06 + 2.466600E+06 -8.400000E-06 + 5.613100E+05 -4.410000E-05 + 1.589900E+05 -1.862000E-04 + 5.186900E+04 -6.783000E-04 + 1.872600E+04 -2.212200E-03 + 7.303600E+03 -6.575200E-03 + 3.029100E+03 -1.793280E-02 + 1.320800E+03 -4.433210E-02 + 6.000300E+02 -9.634780E-02 + 2.819000E+02 -1.696814E-01 + 1.355400E+02 -1.920769E-01 + 6.487000E+01 2.087310E-02 + 3.212900E+01 4.744996E-01 + 1.603700E+01 5.214907E-01 + 7.784900E+00 1.348001E-01 + 3.724700E+00 3.661400E-03 + 1.758300E+00 1.884000E-03 + 5.833100E-01 -4.582000E-04 + 2.785600E-01 2.165000E-04 + 1.182900E-01 -6.780000E-05 + 3 0 21 + 1.647500E+07 4.000000E-07 + 2.466600E+06 3.300000E-06 + 5.613100E+05 1.750000E-05 + 1.589900E+05 7.400000E-05 + 5.186900E+04 2.697000E-04 + 1.872600E+04 8.799000E-04 + 7.303600E+03 2.619800E-03 + 3.029100E+03 7.167100E-03 + 1.320800E+03 1.785610E-02 + 6.000300E+02 3.939600E-02 + 2.819000E+02 7.171020E-02 + 1.355400E+02 8.588770E-02 + 6.487000E+01 -1.038610E-02 + 3.212900E+01 -3.040135E-01 + 1.603700E+01 -4.933178E-01 + 7.784900E+00 1.608900E-02 + 3.724700E+00 7.146686E-01 + 1.758300E+00 4.904795E-01 + 5.833100E-01 3.375730E-02 + 2.785600E-01 -6.792500E-03 + 1.182900E-01 1.691600E-03 + 4 0 1 + 5.833100E-01 1.000000E+00 + 5 0 1 + 2.785600E-01 1.000000E+00 + 6 0 21 + 1.647500E+07 -1.000000E-07 + 2.466600E+06 -1.000000E-06 + 5.613100E+05 -5.400000E-06 + 1.589900E+05 -2.270000E-05 + 5.186900E+04 -8.270000E-05 + 1.872600E+04 -2.694000E-04 + 7.303600E+03 -8.042000E-04 + 3.029100E+03 -2.194900E-03 + 1.320800E+03 -5.493900E-03 + 6.000300E+02 -1.209600E-02 + 2.819000E+02 -2.226230E-02 + 1.355400E+02 -2.660630E-02 + 6.487000E+01 2.758000E-03 + 3.212900E+01 1.016803E-01 + 1.603700E+01 1.704132E-01 + 7.784900E+00 -6.222000E-03 + 3.724700E+00 -3.452570E-01 + 1.758300E+00 -4.234840E-01 + 5.833100E-01 3.354720E-01 + 2.785600E-01 6.779724E-01 + 1.182900E-01 2.523080E-01 + 7 0 1 + 1.182900E-01 1.000000E+00 + 8 0 1 + 0.0442700 1.0000000 + 2 1 16 + 2.660700E+04 6.190000E-05 + 6.298200E+03 5.499000E-04 + 2.045200E+03 3.162000E-03 + 7.821600E+02 1.379790E-02 + 3.316300E+02 4.798120E-02 + 1.511100E+02 1.315710E-01 + 7.239200E+01 2.685861E-01 + 3.586200E+01 3.683473E-01 + 1.813400E+01 2.711363E-01 + 9.043000E+00 7.622220E-02 + 4.450000E+00 4.674900E-03 + 2.166100E+00 1.256500E-03 + 9.962800E-01 -2.357000E-04 + 4.544300E-01 1.332000E-04 + 1.940400E-01 -6.890000E-05 + 7.899700E-02 1.380000E-05 + 3 1 16 + 2.660700E+04 -2.480000E-05 + 6.298200E+03 -2.212000E-04 + 2.045200E+03 -1.273600E-03 + 7.821600E+02 -5.617900E-03 + 3.316300E+02 -1.986000E-02 + 1.511100E+02 -5.655310E-02 + 7.239200E+01 -1.209479E-01 + 3.586200E+01 -1.773098E-01 + 1.813400E+01 -9.214720E-02 + 9.043000E+00 2.187683E-01 + 4.450000E+00 4.854670E-01 + 2.166100E+00 3.721970E-01 + 9.962800E-01 7.769070E-02 + 4.544300E-01 6.541000E-04 + 1.940400E-01 1.799800E-03 + 7.899700E-02 -6.550000E-05 + 4 1 1 + 4.544300E-01 1.000000E+00 + 5 1 16 + 2.660700E+04 6.400000E-06 + 6.298200E+03 5.720000E-05 + 2.045200E+03 3.297000E-04 + 7.821600E+02 1.456200E-03 + 3.316300E+02 5.159100E-03 + 1.511100E+02 1.476170E-02 + 7.239200E+01 3.176940E-02 + 3.586200E+01 4.706800E-02 + 1.813400E+01 2.238710E-02 + 9.043000E+00 -7.202540E-02 + 4.450000E+00 -1.626429E-01 + 2.166100E+00 -1.496503E-01 + 9.962800E-01 1.064517E-01 + 4.544300E-01 4.360266E-01 + 1.940400E-01 4.680523E-01 + 7.899700E-02 1.691513E-01 + 6 1 1 + 1.940400E-01 1.000000E+00 + 7 1 1 + 7.899700E-02 1.000000E+00 + 8 1 1 + 0.0305130 1.0000000 + 3 2 11 + 1.289600E+03 1.190000E-04 + 3.897500E+02 1.155100E-03 + 1.517600E+02 6.764800E-03 + 6.722300E+01 2.730170E-02 + 3.191300E+01 8.092980E-02 + 1.585700E+01 1.794011E-01 + 8.054500E+00 2.840086E-01 + 4.088700E+00 3.266797E-01 + 2.055600E+00 2.584900E-01 + 9.950900E-01 1.128946E-01 + 4.231300E-01 1.726050E-02 + 4 2 1 + 9.950900E-01 1.000000E+00 + 5 2 1 + 4.231300E-01 1.000000E+00 + 6 2 1 + 1.779000E-01 1.000000E+00 + 7 2 1 + 0.0829000 1.0000000 + 4 3 1 + 8.257000E-01 1.000000E+00 + 5 3 1 + 3.407000E-01 1.000000E+00 + 6 3 1 + 0.1748000 1.0000000 + 5 4 1 + 6.491000E-01 1.0000000 + 6 4 1 + 0.3110000 1.0000000 + diff --git a/lib/basis/QZSPDN b/lib/basis/QZSPDN new file mode 100644 index 00000000..eff8c6f8 --- /dev/null +++ b/lib/basis/QZSPDN @@ -0,0 +1,71 @@ +O-HYDROGEN H_1 (QZSPDZDH) BASIS TYPE: 2 +# +10 +1 0 1 +34.63529000 1.00000000 +2 0 1 +23.53722000 1.00000000 +3 0 1 +16.15827000 1.00000000 +4 0 1 +2.88059800 1.00000000 +5 1 1 +33.60130000 1.00000000 +6 1 1 +22.91236000 1.00000000 +7 1 1 +15.80526000 1.00000000 +8 1 1 +6.89794800 1.00000000 +9 2 1 +23.70662000 1.00000000 +10 2 1 +16.73036000 1.00000000 + +O-HYDROGEN H-A_1 (QZSPDZDH) BASIS TYPE: 2 +# +10 +1 0 1 +34.63529000 1.00000000 +2 0 1 +23.53722000 1.00000000 +3 0 1 +16.15827000 1.00000000 +4 0 1 +2.88059800 1.00000000 +5 1 1 +33.60130000 1.00000000 +6 1 1 +22.91236000 1.00000000 +7 1 1 +15.80526000 1.00000000 +8 1 1 +6.89794800 1.00000000 +9 2 1 +23.70662000 1.00000000 +10 2 1 +16.73036000 1.00000000 + +O-HYDROGEN H-B_1 (QZSPDZDH) BASIS TYPE: 2 +# +10 +1 0 1 +34.63529000 1.00000000 +2 0 1 +23.53722000 1.00000000 +3 0 1 +16.15827000 1.00000000 +4 0 1 +2.88059800 1.00000000 +5 1 1 +33.60130000 1.00000000 +6 1 1 +22.91236000 1.00000000 +7 1 1 +15.80526000 1.00000000 +8 1 1 +6.89794800 1.00000000 +9 2 1 +23.70662000 1.00000000 +10 2 1 +16.73036000 1.00000000 diff --git a/src/CI/ArpackInterface.f90 b/src/CI/ArpackInterface.f90 index a7cf55cf..313767ad 100644 --- a/src/CI/ArpackInterface.f90 +++ b/src/CI/ArpackInterface.f90 @@ -22,31 +22,31 @@ !! tomadas de (http://www.netlib.org/blas/) !! - 2011-02-11 : Fernando Posada ( efposadac@unal.edu.co ) !! -# Adapta el modulo para su inclusion en Lowdin -module ArpackInterface_ - implicit none +! module ArpackInterface_ +! implicit none - interface ArpackInterface +! interface ArpackInterface - subroutine dsaupd ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,IPNTR, WORKD, WORKL, LWORKL, INFO ) - character bmat*1, which*2 - integer ido, info, ldv, lworkl, n, ncv, nev - Double precision tol - integer iparam(*), ipntr(*) - Double precision resid(n), v(ldv,ncv), workd(3*n), workl(lworkl) - end subroutine dsaupd +! subroutine dsaupd ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,IPNTR, WORKD, WORKL, LWORKL, INFO ) +! character bmat*1, which*2 +! integer ido, info, ldv, lworkl, n, ncv, nev +! Double precision tol +! integer iparam(*), ipntr(*) +! Double precision resid(n), v(ldv,ncv), workd(3*n), workl(lworkl) +! end subroutine dsaupd - subroutine DSEUPD(RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, & - WORKD, WORKL, LWORKL, INFO ) +! subroutine DSEUPD(RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, & +! WORKD, WORKL, LWORKL, INFO ) - character bmat, howmny, which*2 - logical rvec - integer info, ldz, ldv, lworkl, n, ncv, nev - Double precision sigma, tol - integer iparam(*), ipntr(*) - logical select(ncv) - Double precision d(nev), resid(n), v(ldv,ncv), z(ldz, nev), workd(2*n), workl(lworkl) - end subroutine dseupd +! character bmat, howmny, which*2 +! logical rvec +! integer info, ldz, ldv, lworkl, n, ncv, nev +! Double precision sigma, tol +! integer iparam(*), ipntr(*) +! logical select(ncv) +! Double precision d(nev), resid(n), v(ldv,ncv), z(ldz, nev), workd(2*n), workl(lworkl) +! end subroutine dseupd - end interface +! end interface -end module ArpackInterface_ +! end module ArpackInterface_ diff --git a/src/CI/ConfigurationInteraction.f90 b/src/CI/ConfigurationInteraction.f90 index c1cf1304..10457a9d 100644 --- a/src/CI/ConfigurationInteraction.f90 +++ b/src/CI/ConfigurationInteraction.f90 @@ -39,7 +39,7 @@ module ConfigurationInteraction_ use IndexMap_ use InputCI_ use omp_lib - use ArpackInterface_ + ! use ArpackInterface_ use JadamiluInterface_ implicit none @@ -1943,9 +1943,9 @@ subroutine ConfigurationInteraction_run() select case (trim(String_getUppercase(CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - case ("ARPACK") + ! case ("ARPACK") - write (*,*) "This method was removed" + ! write (*,*) "This method was removed" case ("JADAMILU") diff --git a/src/DFT/Functional.f90 b/src/DFT/Functional.f90 index 270812cf..8b5c7b16 100644 --- a/src/DFT/Functional.f90 +++ b/src/DFT/Functional.f90 @@ -116,6 +116,7 @@ subroutine Functional_constructor(this, speciesID, otherSpeciesID) type(Functional) :: this integer :: speciesID integer :: otherSpeciesID + character(50) :: auxstring this%name="NONE" this%correlationName="NONE" @@ -287,16 +288,20 @@ subroutine Functional_constructor(this, speciesID, otherSpeciesID) !Provisional, only correlation between electron and other species (nuclei) elseif( (this%species1 .eq. "E-" .or. this%species1 .eq. "E-ALPHA" .or. this%species1 .eq. "E-BETA") .and. & (this%species2 .ne. "E-" .and. this%species2 .ne. "E-ALPHA" .and. this%species2 .ne. "E-BETA") ) then + + if (trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring=trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) + else if(trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring=trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) + else + auxstring="NONE" + end if - this%name="correlation:"//trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) - this%correlationName=trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) - this%exactExchangeFraction=1.0_8 !should be irrelevant + this%name="correlation:"//trim(auxstring) + this%correlationName=trim(auxstring) else - this%name="NONE" - this%exactExchangeFraction=1.0_8 - end if diff --git a/src/DFT/GridManager.f90 b/src/DFT/GridManager.f90 index 3c558c8c..16f8c530 100644 --- a/src/DFT/GridManager.f90 +++ b/src/DFT/GridManager.f90 @@ -660,7 +660,7 @@ end subroutine GridManager_getElectronicEnergyAndPotentialAtGrid !!! Felix Moncada, 2017 !< subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid( speciesID, otherSpeciesID, exchangeCorrelationEnergy, & - otherElectronID, otherElectronExchangeCorrelationEnergy) + otherElectronID, otherElectronExchangeCorrelationEnergy) implicit none integer :: speciesID @@ -669,7 +669,7 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid( speciesID, other integer, optional :: otherElectronID real(8), optional :: otherElectronExchangeCorrelationEnergy - character(50) :: nameOfSpecies, otherNameOfSpecies + character(50) :: nameOfSpecies, otherNameOfSpecies, auxstring integer :: gridSize, otherGridSize type(Vector) :: energyDensity type(Vector) :: sigma @@ -680,6 +680,8 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid( speciesID, other nameOfSpecies = MolecularSystem_getNameOfSpecie( speciesID ) otherNameOfSpecies = MolecularSystem_getNameOfSpecie( otherSpeciesID ) + if ( (nameOfSpecies.ne."E-" .and. nameOfSpecies.ne."E-ALPHA") ) return + gridSize = Grid_instance(speciesID)%totalSize otherGridSize = Grid_instance(otherSpeciesID)%totalSize @@ -687,149 +689,153 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid( speciesID, other !Nuclear electron correlation !"E-BETA" is treated at the same time as E-ALPHA - if ( (nameOfSpecies=="E-" .or. nameOfSpecies=="E-ALPHA") ) then - call Vector_constructor(energyDensity, otherGridSize, 0.0_8) - call Vector_constructor(electronicDensityAtOtherGrid, otherGridSize, 0.0_8) - call Vector_constructor(electronicPotentialAtOtherGrid, otherGridSize, 0.0_8) - do dir=1,3 - call Vector_constructor(electronicGradientAtOtherGrid(dir), otherGridSize, 0.0_8) - call Vector_constructor(electronicGradientPotentialAtOtherGrid(dir), otherGridSize, 0.0_8) - end do + call Vector_constructor(energyDensity, otherGridSize, 0.0_8) + call Vector_constructor(electronicDensityAtOtherGrid, otherGridSize, 0.0_8) + call Vector_constructor(electronicPotentialAtOtherGrid, otherGridSize, 0.0_8) + do dir=1,3 + call Vector_constructor(electronicGradientAtOtherGrid(dir), otherGridSize, 0.0_8) + call Vector_constructor(electronicGradientPotentialAtOtherGrid(dir), otherGridSize, 0.0_8) + end do - !!This adds E-BETA density and gradient - call GridManager_getElectronicDensityInOtherGrid(speciesID, otherSpeciesID, & - GridsCommonPoints(speciesID,otherSpeciesID)%totalSize, int(GridsCommonPoints(speciesID,otherSpeciesID)%points%values), & - electronicDensityAtOtherGrid, electronicGradientAtOtherGrid) + !!This adds E-BETA density and gradient + call GridManager_getElectronicDensityInOtherGrid(speciesID, otherSpeciesID, & + GridsCommonPoints(speciesID,otherSpeciesID)%totalSize, int(GridsCommonPoints(speciesID,otherSpeciesID)%points%values), & + electronicDensityAtOtherGrid, electronicGradientAtOtherGrid) + if (trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring=trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) + else if(trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring=trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) + else + auxstring="NONE" + end if !!!These routines return the electronic energy density - select case (trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) ) - - case ("EXPCS-GGA") - call Functional_expCSGGAEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid, electronicGradientAtOtherGrid, & - Grid_instance(otherSpeciesID)%density, Grid_instance(otherSpeciesID)%densityGradient, & - energyDensity, electronicPotentialAtOtherGrid, electronicGradientPotentialAtOtherGrid, & - Grid_instance(otherSpeciesID)%potential, Grid_instance(otherSpeciesID)%gradientPotential ) - - case ("EXPCS-GGA-NOA") - call Functional_expCSGGAEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid, electronicGradientAtOtherGrid, & - Grid_instance(otherSpeciesID)%density, Grid_instance(otherSpeciesID)%densityGradient, & - energyDensity, electronicPotentialAtOtherGrid, electronicGradientPotentialAtOtherGrid, & - Grid_instance(otherSpeciesID)%potential, Grid_instance(otherSpeciesID)%gradientPotential ) - - case ("EPC17-1") - call Functional_EPCEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("EPC17-2") - call Functional_EPCEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("IKN-NSF") - call Functional_IKNEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("MLCS-FIT") - call Functional_MLCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("MLCS-A") - call Functional_MLCSAEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("MLCS-AN") - call Functional_MLCSANEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("CS-MYFIT") - call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("IMAMURA-MYFIT") - call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("MEJIA-MYFIT") - call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("MEJIAA-MYFIT") - call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("EXPCS-A") - call Functional_expCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("EXPCS-NOA") - call Functional_expCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("PSN") - call Functional_PSNEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("PSNAP") - call Functional_PSNAPEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & - electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & - energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) - - case ("NONE") + select case (trim(auxstring) ) + + case ("EXPCS-GGA") + call Functional_expCSGGAEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid, electronicGradientAtOtherGrid, & + Grid_instance(otherSpeciesID)%density, Grid_instance(otherSpeciesID)%densityGradient, & + energyDensity, electronicPotentialAtOtherGrid, electronicGradientPotentialAtOtherGrid, & + Grid_instance(otherSpeciesID)%potential, Grid_instance(otherSpeciesID)%gradientPotential ) + + case ("EXPCS-GGA-NOA") + call Functional_expCSGGAEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid, electronicGradientAtOtherGrid, & + Grid_instance(otherSpeciesID)%density, Grid_instance(otherSpeciesID)%densityGradient, & + energyDensity, electronicPotentialAtOtherGrid, electronicGradientPotentialAtOtherGrid, & + Grid_instance(otherSpeciesID)%potential, Grid_instance(otherSpeciesID)%gradientPotential ) + + case ("EPC17-1") + call Functional_EPCEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("EPC17-2") + call Functional_EPCEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("IKN-NSF") + call Functional_IKNEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("MLCS-FIT") + call Functional_MLCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("MLCS-A") + call Functional_MLCSAEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("MLCS-AN") + call Functional_MLCSANEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("CS-MYFIT") + call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("IMAMURA-MYFIT") + call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("MEJIA-MYFIT") + call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("MEJIAA-MYFIT") + call Functional_myCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("EXPCS-A") + call Functional_expCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("EXPCS-NOA") + call Functional_expCSEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("PSN") + call Functional_PSNEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("PSNAP") + call Functional_PSNAPEvaluate(Functionals(index), MolecularSystem_getMass( otherSpeciesID ), otherGridSize, & + electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values ) + + case ("NONE") + + case default + print *, trim(auxstring) + STOP "The "//otherNameOfSpecies//"electron functional chosen is not implemented" + + end select + + !!Adds the nuclear electron potential to the relevant points in the electronic grid +!!!We partition the energy considering the electronic density + do k=1, GridsCommonPoints(speciesID,otherSpeciesID)%totalSize + !electron index + i=int(GridsCommonPoints(speciesID,otherSpeciesID)%points%values(k,1)) + !nuclear index + j=int(GridsCommonPoints(speciesID,otherSpeciesID)%points%values(k,2)) - case default - print *, trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) - STOP "The nuclear electron functional chosen is not implemented" + exchangeCorrelationEnergy=exchangeCorrelationEnergy+& + energyDensity%values(j)*Grid_instance(speciesID)%density%values(i)*Grid_instance(otherSpeciesID)%points%values(j,4) - end select + Grid_instance(speciesID)%potential%values(i) = Grid_instance(speciesID)%potential%values(i) + electronicPotentialAtOtherGrid%values(j) - !!Adds the nuclear electron potential to the relevant points in the electronic grid -!!!We partition the energy considering the electronic density - do k=1, GridsCommonPoints(speciesID,otherSpeciesID)%totalSize - !electron index - i=int(GridsCommonPoints(speciesID,otherSpeciesID)%points%values(k,1)) - !nuclear index - j=int(GridsCommonPoints(speciesID,otherSpeciesID)%points%values(k,2)) + do dir=1,3 + Grid_instance(speciesID)%gradientPotential(dir)%values(i) = Grid_instance(speciesID)%gradientPotential(dir)%values(i) + electronicGradientPotentialAtOtherGrid(dir)%values(j) + end do - exchangeCorrelationEnergy=exchangeCorrelationEnergy+& - energyDensity%values(j)*Grid_instance(speciesID)%density%values(i)*Grid_instance(otherSpeciesID)%points%values(j,4) + if(nameOfSpecies .eq. "E-ALPHA") then - Grid_instance(speciesID)%potential%values(i) = Grid_instance(speciesID)%potential%values(i) + electronicPotentialAtOtherGrid%values(j) + otherElectronExchangeCorrelationEnergy=otherElectronExchangeCorrelationEnergy+& + energyDensity%values(j)*Grid_instance(otherElectronID)%density%values(i)*Grid_instance(otherSpeciesID)%points%values(j,4) + Grid_instance(otherElectronID)%potential%values(i) = Grid_instance(otherElectronID)%potential%values(i) + electronicPotentialAtOtherGrid%values(j) do dir=1,3 - Grid_instance(speciesID)%gradientPotential(dir)%values(i) = Grid_instance(speciesID)%gradientPotential(dir)%values(i) + electronicGradientPotentialAtOtherGrid(dir)%values(j) + Grid_instance(otherElectronID)%gradientPotential(dir)%values(i) = Grid_instance(otherElectronID)%gradientPotential(dir)%values(i)& + + electronicGradientPotentialAtOtherGrid(dir)%values(j) end do - - if(nameOfSpecies .eq. "E-ALPHA") then - otherElectronExchangeCorrelationEnergy=otherElectronExchangeCorrelationEnergy+& - energyDensity%values(j)*Grid_instance(otherElectronID)%density%values(i)*Grid_instance(otherSpeciesID)%points%values(j,4) - - Grid_instance(otherElectronID)%potential%values(i) = Grid_instance(otherElectronID)%potential%values(i) + electronicPotentialAtOtherGrid%values(j) - do dir=1,3 - Grid_instance(otherElectronID)%gradientPotential(dir)%values(i) = Grid_instance(otherElectronID)%gradientPotential(dir)%values(i)& - + electronicGradientPotentialAtOtherGrid(dir)%values(j) - end do - - end if - - end do + end if - end if + end do call Vector_Destructor(energyDensity) call Vector_destructor(electronicDensityAtOtherGrid) @@ -1059,7 +1065,7 @@ subroutine GridManager_getContactDensity( speciesID, otherSpeciesID, otherElectr integer :: n, nproc real(8) :: contactDensity, overlapDensity - + character(50) :: auxstring type(Vector) :: electronicDensityAtOtherGrid, electronicGradientAtOtherGrid(3), gfactor gridSize =Grid_instance(otherSpeciesID)%totalSize @@ -1084,7 +1090,16 @@ subroutine GridManager_getContactDensity( speciesID, otherSpeciesID, otherElectr ! call Matrix_Constructor( nodeExchangeCorrelationMatrix(n), int(numberOfContractions,8), int(numberOfContractions,8), 0.0_8) ! end do - if(MolecularSystem_getMass(otherSpeciesID) .lt. 2.0 .and. (CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL.eq."expCS-A" .or. CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL.eq."expCS-GGA")) then !positron + if (trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring=trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) + else if(trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring=trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) + else + auxstring="NONE" + end if + + + if(MolecularSystem_getMass(otherSpeciesID) .lt. 2.0 .and. (auxstring.eq."expCS-A" .or. auxstring.eq."expCS-GGA")) then !positron kf=2.2919886876120283056 a0n=0.3647813291441602 @@ -1139,7 +1154,7 @@ subroutine GridManager_getContactDensity( speciesID, otherSpeciesID, otherElectr print *, "Contact density between ", trim(MolecularSystem_getNameOfSpecie(speciesID)),"-", trim(MolecularSystem_getNameOfSpecie(otherSpeciesID)) if(present(otherElectronID) ) print *, "Including contact density between ", trim(MolecularSystem_getNameOfSpecie(otherElectronID)),"-", trim(MolecularSystem_getNameOfSpecie(otherSpeciesID)) - if(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL.eq."expCS-A" .or. CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL.eq."expCS-GGA" ) then + if(auxstring.eq."expCS-A" .or. auxstring.eq."expCS-GGA" ) then print *, "As the integral of rhoA*rhoB(1+g[beta])" print *, "With g[beta] from the expCS-A functional" else diff --git a/src/core/CONTROL.f90 b/src/core/CONTROL.f90 index ad29e403..2d857e48 100644 --- a/src/core/CONTROL.f90 +++ b/src/core/CONTROL.f90 @@ -254,6 +254,7 @@ module CONTROL_ character(50) :: ELECTRON_EXCHANGE_FUNCTIONAL character(50) :: ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL character(50) :: NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL + character(50) :: POSITRON_ELECTRON_CORRELATION_FUNCTIONAL character(50) :: BETA_FUNCTION integer :: GRID_RADIAL_POINTS integer :: GRID_ANGULAR_POINTS @@ -592,6 +593,7 @@ module CONTROL_ character(50) :: LowdinParameters_electronExchangeFunctional character(50) :: LowdinParameters_electronExchangeCorrelationFunctional character(50) :: LowdinParameters_nuclearElectronCorrelationFunctional + character(50) :: LowdinParameters_positronElectronCorrelationFunctional character(50) :: LowdinParameters_betaFunction integer :: LowdinParameters_gridRadialPoints integer :: LowdinParameters_gridAngularPoints @@ -928,6 +930,7 @@ module CONTROL_ LowdinParameters_electronExchangeFunctional,& LowdinParameters_electronExchangeCorrelationFunctional,& LowdinParameters_nuclearElectronCorrelationFunctional,& + LowdinParameters_positronElectronCorrelationFunctional,& LowdinParameters_betaFunction,& LowdinParameters_gridRadialPoints,& LowdinParameters_gridAngularPoints,& @@ -1289,6 +1292,7 @@ subroutine CONTROL_start() LowdinParameters_electronExchangeFunctional = "NONE" LowdinParameters_electronExchangeCorrelationFunctional = "NONE" LowdinParameters_nuclearElectronCorrelationFunctional = "NONE" + LowdinParameters_positronElectronCorrelationFunctional = "NONE" LowdinParameters_betaFunction = "NONE" LowdinParameters_gridRadialPoints=35 LowdinParameters_gridAngularPoints=110 @@ -1625,6 +1629,7 @@ subroutine CONTROL_start() CONTROL_instance%ELECTRON_EXCHANGE_FUNCTIONAL = "NONE" CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL = "NONE" CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL = "NONE" + CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL = "NONE" CONTROL_instance%BETA_FUNCTION = "NONE" CONTROL_instance%GRID_RADIAL_POINTS= 35 CONTROL_instance%GRID_ANGULAR_POINTS= 110 @@ -2018,6 +2023,7 @@ subroutine CONTROL_load(unit) CONTROL_instance%ELECTRON_EXCHANGE_FUNCTIONAL = LowdinParameters_electronExchangeFunctional CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL = LowdinParameters_electronExchangeCorrelationFunctional CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL = LowdinParameters_nuclearElectronCorrelationFunctional + CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL = LowdinParameters_positronElectronCorrelationFunctional CONTROL_instance%BETA_FUNCTION = LowdinParameters_betaFunction CONTROL_instance%GRID_RADIAL_POINTS= LowdinParameters_gridRadialPoints CONTROL_instance%GRID_ANGULAR_POINTS= LowdinParameters_gridAngularPoints @@ -2376,6 +2382,7 @@ subroutine CONTROL_save( unit, lastStep, firstStep ) LowdinParameters_electronExchangeFunctional = CONTROL_instance%ELECTRON_EXCHANGE_FUNCTIONAL LowdinParameters_electronExchangeCorrelationFunctional = CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL LowdinParameters_nuclearElectronCorrelationFunctional = CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL + LowdinParameters_positronElectronCorrelationFunctional = CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL LowdinParameters_betaFunction = CONTROL_instance%BETA_FUNCTION LowdinParameters_gridRadialPoints = CONTROL_instance%GRID_RADIAL_POINTS LowdinParameters_gridAngularPoints = CONTROL_instance%GRID_ANGULAR_POINTS @@ -2702,6 +2709,7 @@ subroutine CONTROL_copy(this, otherThis) otherThis%ELECTRON_EXCHANGE_FUNCTIONAL = this%ELECTRON_EXCHANGE_FUNCTIONAL otherThis%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL = this%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL otherThis%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL = this%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL + otherThis%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL = this%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL otherThis%BETA_FUNCTION = this%BETA_FUNCTION otherThis%GRID_RADIAL_POINTS=this%GRID_RADIAL_POINTS otherThis%GRID_ANGULAR_POINTS=this%GRID_ANGULAR_POINTS @@ -2831,6 +2839,7 @@ subroutine CONTROL_show() end if write (*,"(T10,A)") "ELECTRON-NUCLEAR CORRELATION FUNCTIONAL: "//trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) + write (*,"(T10,A)") "ELECTRON-POSITRON CORRELATION FUNCTIONAL: "//trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) write (*,"(T10,A,I5,A,I5)") "SCF ATOMIC RADIALxANGULAR GRID SIZE:",CONTROL_instance%GRID_RADIAL_POINTS,"x",CONTROL_instance%GRID_ANGULAR_POINTS if( CONTROL_instance%FINAL_GRID_ANGULAR_POINTS*CONTROL_instance%FINAL_GRID_RADIAL_POINTS .gt. & CONTROL_instance%GRID_ANGULAR_POINTS*CONTROL_instance%GRID_RADIAL_POINTS) then diff --git a/src/core/MolecularSystem.f90 b/src/core/MolecularSystem.f90 index e141d5ad..4d4a7b68 100644 --- a/src/core/MolecularSystem.f90 +++ b/src/core/MolecularSystem.f90 @@ -1330,6 +1330,19 @@ function MolecularSystem_getNameOfSpecies(speciesID) result(output) output = MolecularSystem_instance%species(speciesID)%name end function MolecularSystem_getNameOfSpecies + + !> @brief Returns the symbol of a species + !! @author E. F. Posada, 2013 + !! @version 1.0 + function MolecularSystem_getSymbolOfSpecies(speciesID) result(output) + implicit none + + integer :: speciesID + character(30) :: output + + output = MolecularSystem_instance%species(speciesID)%symbol + + end function MolecularSystem_getSymbolOfSpecies !> @brief Returns the name of a species !! @author E. F. Posada, 2013 diff --git a/src/output/OutputBuilder.f90 b/src/output/OutputBuilder.f90 index 8637f48b..b027d18f 100644 --- a/src/output/OutputBuilder.f90 +++ b/src/output/OutputBuilder.f90 @@ -129,7 +129,7 @@ subroutine OutputBuilder_constructor(this, ID, type ,species, state, orbital, di 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) + this%species=MolecularSystem_getNameOfSpecies(1) else allocate(this%fileName(1)) end if @@ -164,7 +164,10 @@ subroutine OutputBuilder_constructor(this, ID, type ,species, state, orbital, di this%auxID=1 !!Check for other outputs of the same type do i=1, this%outputID-1 - if( trim(outputs_instance(i)%type) .eq. trim(this%type) ) this%auxID=this%auxID+1 + if( trim(outputs_instance(i)%type) .eq. trim(this%type) .and. & + trim(outputs_instance(i)%species) .eq. trim(this%species) .and. & + outputs_instance(i)%dimensions .eq. this%dimensions .and. & + outputs_instance(i)%orbital .eq. this%orbital) this%auxID=this%auxID+1 end do end subroutine OutputBuilder_constructor @@ -389,7 +392,7 @@ subroutine OutputBuilder_writeMoldenFile(this) logical :: existFile ! if ( CONTROL_instance%ARE_THERE_DUMMY_ATOMS ) then - ! auxString=MolecularSystem_getNameOfSpecie( 1 ) + ! auxString=MolecularSystem_getNameOfSpecies( 1 ) ! this%fileName=trim(CONTROL_instance%INPUT_FILE)//"mol" ! open(10,file=this%fileName,status='replace',action='write') ! write (10,"(A)") "Hola soy un archivo de molden" @@ -425,7 +428,7 @@ subroutine OutputBuilder_writeMoldenFile(this) write(auxstring,*) state arguments(1) = "OCCUPATIONS"//trim(adjustl(auxstring)) - arguments(2) = MolecularSystem_getNameOfSpecie( l ) + arguments(2) = MolecularSystem_getNameOfSpecies( l ) call Vector_getFromFile(elementsNum=MolecularSystem_getTotalNumberOfContractions(l),& unit=occupationsUnit,& arguments=arguments(1:2),& @@ -461,7 +464,7 @@ subroutine OutputBuilder_writeMoldenFile(this) do i=1, MolecularSystem_getOcupationNumber(l) fractionalOccupations(l,1)%values(i)=1.0_8 * MolecularSystem_getLambda(l) end do - arguments(2) = MolecularSystem_getNameOfSpecie(l) + arguments(2) = MolecularSystem_getNameOfSpecies(l) arguments(1) = "COEFFICIENTS" coefficientsOfcombination(l,1) = & Matrix_getFromFile(unit=wfnUnit, & @@ -488,10 +491,10 @@ subroutine OutputBuilder_writeMoldenFile(this) do l=1,numberOfSpecies if (state .eq. 1) then - auxString=MolecularSystem_getNameOfSpecie( l ) + auxString=MolecularSystem_getNameOfSpecies( l ) else write(auxString, "(I8)") state - auxString=trim(MolecularSystem_getNameOfSpecie( l ))//"-"//trim( adjustl(auxString)) + auxString=trim(MolecularSystem_getNameOfSpecies( l ))//"-"//trim( adjustl(auxString)) end if this%fileName(l)=trim(CONTROL_instance%INPUT_FILE)//trim(auxString)//".molden" @@ -657,7 +660,7 @@ subroutine OutputBuilder_VecGamessFile(this) do l=1,MolecularSystem_getNumberOfQuantumSpecies() - auxString=MolecularSystem_getNameOfSpecie( l ) + auxString=MolecularSystem_getNameOfSpecies( l ) this%fileName(l)=trim(CONTROL_instance%INPUT_FILE)//trim(auxString)//".vec" @@ -665,7 +668,7 @@ subroutine OutputBuilder_VecGamessFile(this) specieID = int( MolecularSystem_getSpecieID(nameOfSpecie = trim(auxString)) ) numberOfContractions = MolecularSystem_getTotalNumberOfContractions(specieID) - arguments(2) = MolecularSystem_getNameOfSpecie(specieID) + arguments(2) = MolecularSystem_getNameOfSpecies(specieID) arguments(1) = "COEFFICIENTS" coefficientsOfcombination = & @@ -742,36 +745,21 @@ end subroutine OutputBuilder_VecGamessFile subroutine OutputBuilder_casinoFile(this) implicit none type(OutputBuilder) :: this - type(MolecularSystem) :: MolecularSystemInstance integer :: i integer :: j - integer :: k integer :: l integer :: g, h, m integer :: specieID - logical :: wasPress - character(10) :: auxString - character(10) :: symbol - real(8) :: origin(3) - real(8), allocatable :: charges(:) - type(Matrix) :: localizationOfCenters - type(Matrix) :: auxMatrix type(Matrix) :: coefficientsOfcombination real(8), allocatable :: superMatrix(:,:) - character(10),allocatable :: labels(:) integer :: wfnUnit character(50) :: wfnFile integer :: numberOfContractions, superSize integer :: numberOfContractionsA, numberOfContractionsB integer :: numberOfShellsA, numberOfShellsB, totalShells character(50) :: arguments(20) - character(19) , allocatable :: labelsOfContractions(:) - integer :: counter, auxcounter - character(6) :: nickname - character(2) :: space integer :: i0, j0, maxl, shellCode - integer :: totalNumberOfParticles, n real(8) :: puntualInteractionEnergy wfnFile = "lowdin.wfn" @@ -1038,7 +1026,7 @@ subroutine OutputBuilder_casinoFile(this) specieID = l numberOfContractions = MolecularSystem_getTotalNumberOfContractions(specieID) - arguments(2) = MolecularSystem_getNameOfSpecie(specieID) + arguments(2) = MolecularSystem_getNameOfSpecies(specieID) arguments(1) = "COEFFICIENTS" coefficientsOfcombination = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2)) @@ -1217,7 +1205,7 @@ subroutine OutputBuilder_writeEigenvalues(this) wfnFile = "lowdin.wfn" wfnUnit = 20 - ! auxString=MolecularSystem_getNameOfSpecie( 1 ) + ! auxString=MolecularSystem_getNameOfSpecies( 1 ) ! this%fileName=trim(CONTROL_instance%INPUT_FILE)//".eigen" ! open(129,file=this%fileName,status='replace',action='write') ! close(129) @@ -1237,14 +1225,14 @@ subroutine OutputBuilder_writeEigenvalues(this) totalNumberOfParticles = 0 - auxString=MolecularSystem_getNameOfSpecie( l ) + auxString=MolecularSystem_getNameOfSpecies( l ) specieID = MolecularSystem_getSpecieID(auxString) this%fileName(l)=trim(CONTROL_instance%INPUT_FILE)//trim(auxString)//".eigen" open(129,file=this%fileName(l),status='replace',action='write') specieID = int( MolecularSystem_getSpecieID(nameOfSpecie = trim(auxString)) ) numberOfContractions = MolecularSystem_getTotalNumberOfContractions(specieID) - arguments(2) = MolecularSystem_getNameOfSpecie(specieID) + arguments(2) = MolecularSystem_getNameOfSpecies(specieID) arguments(1) = "ORBITALS" call Vector_getFromFile( elementsNum = numberOfContractions, & @@ -1252,7 +1240,7 @@ subroutine OutputBuilder_writeEigenvalues(this) output = energyOfMolecularOrbital ) do j=1,size(energyOfMolecularOrbital%values) - write (129,"(F15.12)") ,energyOfMolecularOrbital%values(j) + write (129,"(F15.12)") energyOfMolecularOrbital%values(j) end do close(129) end do @@ -1320,7 +1308,7 @@ subroutine OutputBuilder_writeFchkFile(this) ! open(unit = occupationsUnit, file=trim(occupationsFile), status="old", form="formatted") ! do l=1,numberOfSpecies ! arguments(1) = "OCCUPATIONS" - ! arguments(2) = MolecularSystem_getNameOfSpecie( l ) + ! arguments(2) = MolecularSystem_getNameOfSpecies( l ) ! fractionalOccupations(l)= Matrix_getFromFile(unit=occupationsUnit,& ! rows=int(MolecularSystem_getTotalNumberOfContractions(l),4),& ! columns=int(numberOfStates,4),& @@ -1345,7 +1333,7 @@ subroutine OutputBuilder_writeFchkFile(this) ! do state=1,numberOfStates do l=1,numberOfSpecies - nameOfSpecies=MolecularSystem_getNameOfSpecie(l) + nameOfSpecies=MolecularSystem_getNameOfSpecies(l) particlesPerOrbital=MolecularSystem_getLambda(l) ! if (state .eq. 1) then ! auxString=nameOfSpecies @@ -1361,7 +1349,7 @@ subroutine OutputBuilder_writeFchkFile(this) numberOfShells=MolecularSystem_getNumberOfContractions(l) numberOfContractions=MolecularSystem_getTotalNumberOfContractions(l) - arguments(2) = MolecularSystem_getNameOfSpecie(l) + arguments(2) = MolecularSystem_getNameOfSpecies(l) arguments(1) = "COEFFICIENTS" coefficientsOfcombination = & Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), & @@ -1741,7 +1729,7 @@ subroutine OutputBuilder_generateAIMFiles (this) close(35) do l=1,MolecularSystem_getNumberOfQuantumSpecies() - auxString=MolecularSystem_getNameOfSpecie( l ) + auxString=MolecularSystem_getNameOfSpecies( l ) moldenFileName=trim(CONTROL_instance%INPUT_FILE)//trim(auxString)//".molden" call Molden2AIM(moldenFileName, totalEnergy, virial) !! Just for printing information @@ -1760,7 +1748,7 @@ subroutine OutputBuilder_generateExtendedWfnFile (this) character(50) :: auxString do l=1,MolecularSystem_getNumberOfQuantumSpecies() - auxString=MolecularSystem_getNameOfSpecie( l ) + auxString=MolecularSystem_getNameOfSpecies( l ) end do end subroutine OutputBuilder_generateExtendedWfnFile @@ -1778,10 +1766,11 @@ subroutine OutputBuilder_get3DPlot(this) type(vector) :: step2 real(8) :: maxValue, maxValue2 real(8) :: minValue, minValue2 - Type(Vector) :: val, val2 + real(8) :: plotDistance1, plotDistance2 + Type(Vector) :: val Type(Matrix) :: coordinate - character(50) :: title, title2 + character(50) :: title character(50) :: x_title character(50) :: y_title character(50) :: z_title @@ -1789,10 +1778,12 @@ subroutine OutputBuilder_get3DPlot(this) call Vector_Constructor(step1, 3) call Vector_Constructor(step2, 3) speciesID = MolecularSystem_getSpecieIDFromSymbol( trim(this%species) ) - nameOfSpecies=MolecularSystem_getNameOfSpecie(speciesID) + nameOfSpecies=MolecularSystem_getNameOfSpecies(speciesID) this%fileName2="" numberOfSteps= CONTROL_instance%NUMBER_OF_POINTS_PER_DIMENSION + plotDistance1=sqrt(sum((this%point2%values(:)-this%point1%values(:))**2)) + plotDistance2=sqrt(sum((this%point3%values(:)-this%point1%values(:))**2)) step1%values(:)=(this%point2%values(:)-this%point1%values(:))/numberOfSteps step2%values(:)=(this%point3%values(:)-this%point1%values(:))/numberOfSteps outputID=String_convertIntegerToString(this%outputID) @@ -1806,20 +1797,24 @@ subroutine OutputBuilder_get3DPlot(this) case ( "ORBITALPLOT") orbitalNum=String_convertIntegerToString(this%orbital) - this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".3D.orb"//trim(orbitalNum) + if( this%auxID .eq. 1) then + this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".3D.orb"//trim(orbitalNum) + else + this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".3D-"//trim(auxID)//".orb"//trim(orbitalNum) + end if open(10,file=this%fileName(1),status='replace',action='write') write (10,"(A10,A20,A20,A20)") "#", "X","Y","OrbitalValue" title=trim(this%species)//" Orbital Number: "//trim(orbitalNum) - case ( "FUKUIPLOT") - this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".3D.fkpos" - this%fileName2=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".3D.fkneg" - open(10,file=this%fileName(1),status='replace',action='write') - write (10,"(A10,A20,A20,A20)") "#", "X","Y","PositiveFukuiValue" - title=trim(this%species)//" Positive Fukui" - open(11,file=this%fileName2,status='replace',action='write') - write (11,"(A10,A20,A20,A20)") "#", "X","Y","NegativeFukuiValue" - title2=trim(this%species)//" Negative Fukui" + ! case ( "FUKUIPLOT") + ! this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".3D.fkpos" + ! this%fileName2=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".3D.fkneg" + ! open(10,file=this%fileName(1),status='replace',action='write') + ! write (10,"(A10,A20,A20,A20)") "#", "X","Y","PositiveFukuiValue" + ! title=trim(this%species)//" Positive Fukui" + ! open(11,file=this%fileName2,status='replace',action='write') + ! write (11,"(A10,A20,A20,A20)") "#", "X","Y","NegativeFukuiValue" + ! title2=trim(this%species)//" Negative Fukui" case default call OutputBuilder_exception(ERROR, "The output plot type you requested has not been implemented yet", "OutputBuilder_get3DPlot" ) @@ -1842,7 +1837,7 @@ subroutine OutputBuilder_get3DPlot(this) select case( this%type ) case ( "ORBITALPLOT") call CalculateWaveFunction_getOrbitalValueAt( speciesID, this%orbital, coordinate, val ) - case ( "FUKUIPLOT") + ! case ( "FUKUIPLOT") !! val=CalculateProperties_getFukuiAt( this%species, "positive", coordinate ) !! val2=CalculateProperties_getFukuiAt( this%species, "negative", coordinate ) case default @@ -1851,28 +1846,38 @@ subroutine OutputBuilder_get3DPlot(this) n=0 do i=0,numberOfSteps write (10,*) "" - if (this%type .eq. "FUKUIPLOT") write(11,*) "" + ! if (this%type .eq. "FUKUIPLOT") write(11,*) "" do j=0,numberOfSteps n=n+1 - write (10,"(T10,F20.8,F20.8,ES20.8)") i*Vector_norm(step1),j*Vector_norm(step2),val%values(n) + if(abs(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) - if (this%type .eq. "FUKUIPLOT" ) then - 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 + ! if (this%type .eq. "FUKUIPLOT" ) then + ! 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 ! print *, coordinate, val end do end do + !!large orbital values lead to bad looking plots + if(maxValue .gt. 1.0) maxValue=1.0 + if(minValue .lt. -1.0) minValue=-1.0 + call OutputBuilder_make3DGraph( this%fileName(1), title, x_title, y_title, z_title, minValue, maxValue) close(10) - if (this%type .eq. "FUKUIPLOT" ) then - call OutputBuilder_make3DGraph( this%fileName2, title2, x_title, y_title, z_title, minValue2, maxValue2) - close(11) - end if + ! if (this%type .eq. "FUKUIPLOT" ) then + ! call OutputBuilder_make3DGraph( this%fileName2, title2, x_title, y_title, z_title, minValue2, maxValue2) + ! close(11) + ! end if end subroutine OutputBuilder_get3DPlot @@ -1886,21 +1891,23 @@ subroutine OutputBuilder_get2DPlot(this) character(50) :: nameOfSpecies integer :: numberOfSteps - type(vector) :: step - type(vector) :: val, val2 + type(vector) :: step1 + type(vector) :: val Type(Matrix) :: coordinate + real(8) :: plotDistance1 character(50) :: title character(50) :: x_title character(50) :: y_title - call Vector_Constructor(step, 3) + call Vector_Constructor(step1, 3) speciesID = MolecularSystem_getSpecieIDFromSymbol( trim(this%species) ) - nameOfSpecies=MolecularSystem_getNameOfSpecie(speciesID) + nameOfSpecies=MolecularSystem_getNameOfSpecies(speciesID) this%fileName2="" numberOfSteps= CONTROL_instance%NUMBER_OF_POINTS_PER_DIMENSION - step%values(:)=(this%point2%values(:)-this%point1%values(:))/numberOfSteps + plotDistance1=sqrt(sum((this%point2%values(:)-this%point1%values(:))**2)) + step1%values(:)=(this%point2%values(:)-this%point1%values(:))/numberOfSteps outputID=String_convertIntegerToString(this%outputID) auxID=String_convertIntegerToString(this%auxID) @@ -1909,23 +1916,27 @@ subroutine OutputBuilder_get2DPlot(this) case ( "ORBITALPLOT") orbitalNum=String_convertIntegerToString(this%orbital) - this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".2D.orb"//trim(orbitalNum) + if( this%auxID .eq. 1) then + this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".2D.orb"//trim(orbitalNum) + else + this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//trim(nameOfSpecies)//".2D-"//trim(auxID)//".orb"//trim(orbitalNum) + end if open(10,file=this%fileName(1),status='replace',action='write') write (10,"(A10,A20,A20)") "#", "X","OrbitalValue" title=trim(this%species)//" Orbital Number "//trim(orbitalNum) y_title="orbitalValue/a.u.^{-3/2}" - case ( "FUKUIPLOT") - this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".2D.fkpos" - this%fileName2=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".2D.fkneg" + ! case ( "FUKUIPLOT") + ! this%fileName(1)=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".2D.fkpos" + ! this%fileName2=trim(CONTROL_instance%INPUT_FILE)//"out"//trim(outputID)//"."//trim(this%species)//".2D.fkneg" - open(10,file=this%fileName(1),status='replace',action='write') - write (10,"(A10,A20,A20)") "#","X","PositiveFukuiValue" - title=trim(this%species)//" positive fukui" - y_title="density/a.u.^{-3}" + ! open(10,file=this%fileName(1),status='replace',action='write') + ! write (10,"(A10,A20,A20)") "#","X","PositiveFukuiValue" + ! title=trim(this%species)//" positive fukui" + ! y_title="density/a.u.^{-3}" - open(11,file=this%fileName2,status='replace',action='write') - write (11,"(A10,A20,A20)") "#","X","NegativeFukuiValue" + ! open(11,file=this%fileName2,status='replace',action='write') + ! write (11,"(A10,A20,A20)") "#","X","NegativeFukuiValue" case default call OutputBuilder_exception(ERROR, "The output plot type you requested has not been implemented yet", "OutputBuilder_get3DPlot" ) @@ -1934,13 +1945,13 @@ subroutine OutputBuilder_get2DPlot(this) call Matrix_constructor(coordinate,int((numberOfSteps+1),8),int(3,8),0.0_8) do i=0,numberOfSteps - coordinate%values(i+1,:)=i*step%values(:)+this%point1%values(:) + coordinate%values(i+1,:)=i*step1%values(:)+this%point1%values(:) end do select case( this%type ) case ( "ORBITALPLOT") call CalculateWaveFunction_getOrbitalValueAt( speciesID, this%orbital, coordinate, val ) - case ( "FUKUIPLOT") + ! case ( "FUKUIPLOT") !! val=CalculateProperties_getFukuiAt( this%species, "positive", coordinate ) !! val2=CalculateProperties_getFukuiAt( this%species, "negative", coordinate ) case default @@ -1949,8 +1960,12 @@ subroutine OutputBuilder_get2DPlot(this) n=0 do i=0,numberOfSteps n=n+1 - 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) + if(abs(val%values(n))>1.0E-99_8) then + write (10,"(T10,F20.8,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1),val%values(n) + else + write (10,"(T10,F20.8,F20.8,E20.8)") -plotDistance1*0.5+i*Vector_norm(step1),0.0 + end if + ! if (this%type .eq. "FUKUIPLOT") write (11,"(T10,F20.8,ES20.8)") i*Vector_norm(step),val2%values(n) end do close(10) @@ -1961,7 +1976,7 @@ subroutine OutputBuilder_get2DPlot(this) !! title=trim(this%species)//" negative fukui" !! call OutputBuilder_make2DGraph( this%fileName2, title, x_title, y_title) !! end if - call Vector_Destructor ( step) + call Vector_Destructor ( step1) end subroutine OutputBuilder_get2DPlot @@ -1983,7 +1998,7 @@ subroutine OutputBuilder_getDensityCube(this ) integer :: numberOfOrbitals, numberOfSpecies type(matrix) :: densityMatrix - character(100) :: arguments(2), wfnFile, occupationsFile, auxstring, nameOfSpecies + character(100) :: arguments(2), wfnFile, occupationsFile, auxstring, nameOfSpecies, symbolOfSpecies logical :: existFile !Writes Gaussian Cube @@ -1992,8 +2007,9 @@ subroutine OutputBuilder_getDensityCube(this ) l=0 do speciesID=1, numberOfSpecies - nameOfSpecies=MolecularSystem_getNameOfSpecie(speciesID) - if(trim(this%species) .eq. trim(nameOfSpecies) .or. trim(this%species) .eq. "ALL" ) then + nameOfSpecies=MolecularSystem_getNameOfSpecies(speciesID) + symbolOfSpecies=MolecularSystem_getSymbolOfSpecies(speciesID) + if(trim(this%species) .eq. trim(symbolOfSpecies) .or. trim(this%species) .eq. "ALL" ) then l=l+1 numberOfOrbitals=MolecularSystem_getTotalNumberOfContractions(speciesID) @@ -2107,7 +2123,7 @@ subroutine OutputBuilder_getDensityPlot(this) Type(Vector) :: val Type(Matrix) :: coordinate - character(100) :: arguments(2), wfnFile, occupationsFile, auxstring, nameOfSpecies + character(100) :: arguments(2), wfnFile, occupationsFile, auxstring, nameOfSpecies, symbolOfSpecies character(50) :: title, x_title, y_title, z_title logical :: existFile @@ -2120,8 +2136,9 @@ subroutine OutputBuilder_getDensityPlot(this) l=0 do speciesID=1, numberOfSpecies - nameOfSpecies=MolecularSystem_getNameOfSpecie(speciesID) - if(trim(this%species) .eq. trim(nameOfSpecies) .or. trim(this%species) .eq. "ALL" ) then + nameOfSpecies=MolecularSystem_getNameOfSpecies(speciesID) + symbolOfSpecies=MolecularSystem_getSymbolOfSpecies(speciesID) + if(trim(this%species) .eq. trim(symbolOfSpecies) .or. trim(this%species) .eq. "ALL" ) then l=l+1 numberOfOrbitals=MolecularSystem_getTotalNumberOfContractions(speciesID) @@ -2222,7 +2239,7 @@ subroutine OutputBuilder_getDensityPlot(this) end do !!large density values lead to bad looking plots - if(maxValue .gt. 0.5) maxValue=0.5 + if(maxValue .gt. 1.0) maxValue=1.0 call OutputBuilder_make3DGraph( this%fileName(l), title, x_title, y_title, z_title, 0.0_8, maxValue) close(10) @@ -2348,12 +2365,6 @@ subroutine OutputBuilder_make3DGraph(fileName, title, x_title, y_title, z_title, character(*) :: z_title real(8) :: minValue real(8) :: maxValue - real(8) :: maxMinDiff - - integer :: levels - - maxMinDiff=maxValue-minValue - levels=10 open ( 100,FILE=trim(fileName)//'.gnp', STATUS='REPLACE',ACTION='WRITE') write (100,"(A)") 'set term post eps enh color "Helvetica" 16 size 7cm,5cm' @@ -2364,12 +2375,15 @@ subroutine OutputBuilder_make3DGraph(fileName, title, x_title, y_title, z_title, write (100,"(A)") 'splot "'//trim(fileName)//'" u 1:2:3' write (100,"(A)") 'unset table' - write (100,"(A,I5)") 'levels=', levels - write (100,"(A,I5)") 'numColors=', 5 + if(minValue.lt.0 .and. maxValue.gt.0) then + write (100,"(A,I5)") 'levels=', 11 + else + write (100,"(A,I5)") 'levels=', 10 + end if + write (100,"(A,E20.8)") 'maxValue=', maxValue write (100,"(A,E20.8)") 'minValue=', minValue write (100,"(A)") 'step=(maxValue-minValue)/levels' - write (100,"(A)") 'colorStep=(maxValue-minValue)/numColors' write (100,"(A)") 'set contour base' write (100,"(A)") 'set cntrparam level incremental minValue, step , maxValue' @@ -2384,12 +2398,22 @@ subroutine OutputBuilder_make3DGraph(fileName, title, x_title, y_title, z_title, write (100,"(A)") 'set cbrange [minValue:maxValue]' write (100,"(A)") 'set palette maxcolors levels' - write (100,"(A)") 'set cbtics step' + if(minValue.lt.0 .and. maxValue.gt.0) then + write (100,"(A)") 'set cbtics (minValue, 0.0, maxValue)' + else + write (100,"(A)") 'set cbtics step' + end if write (100,"(A)") 'set format cb "%3.1E"' - write (100,"(A)") 'set palette defined (minValue "white",minValue+colorStep "blue",minValue+colorStep*2 "green",minValue+colorStep*3 "yellow",minValue+colorStep*4 "orange",maxValue "red")' - write (100,"(A)") 'set grid front' + if(minValue.lt.0 .and. maxValue.gt.0) then + write (100,"(A)") 'set palette defined (minValue "blue", 0.0 "white", maxValue "red")' + else if(minValue.ge.0) then + write (100,"(A)") 'set palette defined (minValue "white", maxValue "red")' + else + write (100,"(A)") 'set palette defined (minValue "blue", maxValue "white")' + end if + write (100,"(A)") 'set grid front' write (100,"(A)") 'set format x "%.0f"' write (100,"(A)") 'set format y "%.0f"' @@ -2397,71 +2421,7 @@ subroutine OutputBuilder_make3DGraph(fileName, title, x_title, y_title, z_title, write (100,"(A)") 'set ylabel "Y (a.u.)"' write (100,"(A)") 'plot "'//trim(fileName)//'.table" with image, "'//trim(fileName)//'.cont" w l lt -1 lw 1.5' - - - ! write (100,"(A)") 'set output "'//trim(fileName)//'.eps"' - ! write (100,"(A)") 'set xlabel "'//trim(x_title)//'"' - ! write (100,"(A)") 'set ylabel "'//trim(y_title)//'"' - ! write (100,"(A)") 'set zlabel "'//trim(z_title)//'"' - ! write (100,"(A)") 'set pm3d ' - - ! if (minValue < CONTROL_instance%DOUBLE_ZERO_THRESHOLD .and. maxValue > -CONTROL_instance%DOUBLE_ZERO_THRESHOLD) then - ! write (100,"(A)") 'set palette model RGB defined ('//String_convertRealToString(minValue)//& - ! ' "violet", '//String_convertRealToString(minValue/5)//& - ! ' "blue", '//String_convertRealToString(minValue/25)//& - ! ' "green", '//String_convertRealToString(0.0_8) // & - ! ' "white", '//String_convertRealToString(maxValue/25)//& - ! ' "yellow", '//String_convertRealToString(maxValue/5)//& - ! ' "orange", '//String_convertRealToString(maxValue)//' "red") ' - ! end if - - ! if (minValue > -CONTROL_instance%DOUBLE_ZERO_THRESHOLD) then - ! write (100,"(A)") 'set palette model RGB defined ('//String_convertRealToString(maxValue/15625) // & - ! ' "white", '//String_convertRealToString(maxValue/3125)//& - ! ' "violet", '//String_convertRealToString(maxValue/625)//& - ! ' "blue", '//String_convertRealToString(maxValue/125)//& - ! ' "green", '//String_convertRealToString(maxValue/25)//& - ! ' "yellow", '//String_convertRealToString(maxValue/5)//& - ! ' "orange", '//String_convertRealToString(maxValue)//' "red") ' - ! end if - - ! if (maxValue < CONTROL_instance%DOUBLE_ZERO_THRESHOLD) then - ! write (100,"(A)") 'set palette model RGB defined ('//String_convertRealToString(minValue/15625) // & - ! ' "white", '//String_convertRealToString(minValue/3125)//& - ! ' "violet", '//String_convertRealToString(minValue/625)//& - ! ' "blue", '//String_convertRealToString(minValue/125)//& - ! ' "green", '//String_convertRealToString(minValue/25)//& - ! ' "yellow", '//String_convertRealToString(minValue/5)//& - ! ' "orange", '//String_convertRealToString(minValue)//' "red") ' - ! end if - - ! write (100,"(A)") 'set format "%.0f"' - ! write (100,"(A)") 'set format cb "%.1f"' - ! write (100,"(A)") 'set cbtics '//String_convertRealToString(maxValue/5) - ! write (100,"(A)") 'set xyplane at 0' - ! write (100,"(A)") 'set surface' - ! write (100,"(A)") 'set border 4095' - ! write (100,"(A)") 'set cntrparam cubicspline' - ! write (100,"(A)") 'set cntrparam points 20' - ! write (100,"(A)") 'set cntrparam levels 10' - ! write (100,"(A)") 'set rmargin -1' - ! write (100,"(A)") 'set lmargin -1' - ! write (100,"(A)") 'set tmargin -1' - ! write (100,"(A)") 'set bmargin -1' - ! write (100,"(A)") 'unset ztics' - ! write (100,"(A)") 'set multiplot title "'//trim(title)//'" layout 1,2' - - ! write (100,"(A)") 'set colorbox vertical' - ! write (100,"(A)") 'set colorbox user origin 0.48,0.25 size 0.04,0.5' - ! write (100,"(A)") 'set view 50,160' - ! write (100,"(A)") 'splot "'//trim(fileName)//'" u 1:2:3 notitle w pm3d' - - ! write (100,"(A)") 'unset colorbox' - ! write (100,"(A)") 'set view 0,0' - ! write (100,"(A)") 'splot "'//trim(fileName)//'" u 1:2:3 notitle w pm3d' - - ! write (100,"(A)") 'unset multiplot' - + close(100) ! status= system("gnuplot "//trim(fileName)//".gnp") diff --git a/test/CO.CISD-dipole.lowdin b/test/CO.CISD-dipole.lowdin new file mode 100644 index 00000000..889f98a1 --- /dev/null +++ b/test/CO.CISD-dipole.lowdin @@ -0,0 +1,30 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Test CISD +% Electronic Hartree Fock +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SYSTEM_DESCRIPTION='CO' + +GEOMETRY +e-(C) 6-31G.D 0.000 0.000 -0.914059 +e-(O) 6-31G.D 0.000 0.000 0.213941 +C dirac 0.000 0.000 -0.914059 +O dirac 0.000 0.000 0.213941 +END GEOMETRY + +TASKS + method = "UHF" + configurationInteractionLevel ="CISD" +END TASKS + +CONTROL + readCoefficients=.F. + totalenergytolerance=1E-10 + CIdiagonalizationMethod="JADAMILU" +END CONTROL +INPUT_CI + species="E-ALPHA" core = 2 active = 20 + species="E-BETA" core = 2 active = 20 +END INPUT_CI + + diff --git a/test/PsOH.HF.densPlot.py b/test/CO.CISD-dipole.py similarity index 52% rename from test/PsOH.HF.densPlot.py rename to test/CO.CISD-dipole.py index b76b7835..be65eb59 100644 --- a/test/PsOH.HF.densPlot.py +++ b/test/CO.CISD-dipole.py @@ -12,14 +12,13 @@ testName = sys.argv[0][:-3] inputName = testName + ".lowdin" outputName = testName + ".out" -plot1Name = testName + ".E-.3D.dens" -plot2Name = testName + ".POSITRON.3D.dens" # Reference values and tolerance refValues = { -"HF energy" : [-75.587683637788,1E-8], -"Num e- in plot" : [10.45449875,1E-1], -"Num e+ in plot" : [0.99246436,1E-2], +"HF energy" : [-112.737336707933,1E-8], +"CISD energy" : [-112.957030012578,1E-6], +"HF z-dipole" : [-0.33130728,1E-3], +"CISD z-dipole" : [0.13640258,1E-3], } testValues = dict(refValues) #copy @@ -38,39 +37,29 @@ outputRead = output.readlines() # Values +dipoleflag=False +ciflag=False +hfflag=True for i in range(0,len(outputRead)): line = outputRead[i] if "TOTAL ENERGY =" in line: testValues["HF energy"] = float(line.split()[3]) + if "STATE: 1 ENERGY =" in line: + testValues["CISD energy"] = float(line.split()[4]) + if "DIPOLE: (DEBYE)" in line: + dipoleflag=True + if "Total Dipole:" in line and dipoleflag and hfflag: + testValues["HF z-dipole"] = float(line.split()[4]) + dipoleflag=False + hfflag=False + if "We are calculating properties for E-ALPHA in the CI ground state" in line: + ciflag=True + if "Total Dipole:" in line and dipoleflag and ciflag: + testValues["CISD z-dipole"] = float(line.split()[4]) + dipoleflag=False output.close() -plot1 = open(plot1Name, "r") -plot1Read = plot1.readlines() -sumE=0 -for i in range(0,len(plot1Read)): - line = plot1Read[i] - if i > 1: - values = line.split() - if i == 3: x1=float(values[1]) - if i == 4: x2=float(values[1]) - if len(values) > 1: sumE+=(float(values[0])**2+float(values[1])**2)**(1.0/2.0)*float(values[2]) -testValues["Num e- in plot"]=2.0*sumE*(x2-x1)**2 -plot1.close() - -plot2 = open(plot2Name, "r") -plot2Read = plot2.readlines() -sumP=0 -for i in range(0,len(plot2Read)): - line = plot2Read[i] - if i > 1: - values = line.split() - if i == 3: x1=float(values[1]) - if i == 4: x2=float(values[1]) - if len(values) > 1: sumP+=(float(values[0])**2+float(values[1])**2)**(1.0/2.0)*float(values[2]) -testValues["Num e+ in plot"]=2.0*sumP*(x2-x1)**2 -plot2.close() - passTest = True for value in refValues: diff --git a/test/Gly.e+-molden.lowdin b/test/Gly.e+-molden.lowdin new file mode 100644 index 00000000..7f91647a --- /dev/null +++ b/test/Gly.e+-molden.lowdin @@ -0,0 +1,37 @@ +SYSTEM_DESCRIPTION='glicina con postiron' + +GEOMETRY +e-(N) 6-311G 0.92914 0.65159 0.25796 +e-(C) 6-311G 0.08973 -0.44202 0.82116 +e-(C) 6-311G 0.09416 -0.39610 2.35537 +e-(O) 6-311G 0.79215 0.47431 2.87347 +e-(O) 6-311G -0.60778 -1.24843 2.89495 +e-(H) 6-311G 1.89448 0.53627 0.52036 +e-(H) 6-311G 0.48806 -1.38246 0.47040 +e-(H) 6-311G -0.91478 -0.32058 0.44545 +e+ PSX-TZ -0.60778 -1.24843 2.89495 +N dirac 0.92914 0.65159 0.25796 +C dirac 0.08973 -0.44202 0.82116 +C dirac 0.09416 -0.39610 2.35537 +O dirac 0.79215 0.47431 2.87347 +O dirac -0.60778 -1.24843 2.89495 +H dirac 1.89448 0.53627 0.52036 +H dirac 0.48806 -1.38246 0.47040 +H dirac -0.91478 -0.32058 0.44545 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=F +END CONTROL + +OUTPUTS + moldenFile +END OUTPUTS + + + + diff --git a/test/PsF.HF.densPlot.py b/test/Gly.e+-molden.py similarity index 53% rename from test/PsF.HF.densPlot.py rename to test/Gly.e+-molden.py index 57b0c83b..eba52463 100644 --- a/test/PsF.HF.densPlot.py +++ b/test/Gly.e+-molden.py @@ -12,14 +12,14 @@ testName = sys.argv[0][:-3] inputName = testName + ".lowdin" outputName = testName + ".out" -plot1Name = testName + ".E-.2D.dens" -plot2Name = testName + ".POSITRON.2D.dens" +molden1Name = testName + ".E-.molden" +molden2Name = testName + ".POSITRON.molden" # Reference values and tolerance refValues = { -"HF energy" : [-99.635031860198,1E-8], -"Num e- in plot" : [9.99830528,1E-3], -"Num e+ in plot" : [0.99999922,1E-5], +"HF energy" : [-281.379249259318,1E-8], +"e+HOMO" : [-3.48415E-03,1E-4], +"e-HOMO" : [-5.40431E-01,1E-4], } testValues = dict(refValues) #copy @@ -45,31 +45,33 @@ output.close() -plot1 = open(plot1Name, "r") -plot1Read = plot1.readlines() -sumE=0 -for i in range(0,len(plot1Read)): - line = plot1Read[i] - if i > 1: - values = line.split() - if i == 3: x1=float(values[0]) - if i == 4: x2=float(values[0]) - if len(values) > 1: sumE+=float(values[0])**2*float(values[1]) -testValues["Num e- in plot"]=2.0*3.14159265359*sumE*(x2-x1) -plot1.close() +molden1 = open(molden1Name, "r") +molden1Read = molden1.readlines() +v=0 +eigenv=[] +for i in range(0,len(molden1Read)): + line = molden1Read[i] + if "Ene=" in line: + eigenv.append(float(line.split()[1])) + v+=1 + if "Occup=" in line and "0.0" in line : + testValues["e-HOMO"] = eigenv[v-2] + break +molden1.close() -plot2 = open(plot2Name, "r") -plot2Read = plot2.readlines() -sumP=0 -for i in range(0,len(plot2Read)): - line = plot2Read[i] - if i > 1: - values = line.split() - if i == 3: x1=float(values[0]) - if i == 4: x2=float(values[0]) - if len(values) > 1: sumP+=float(values[0])**2*float(values[1]) -testValues["Num e+ in plot"]=2.0*3.14159265359*sumP*(x2-x1) -plot2.close() +molden2 = open(molden2Name, "r") +molden2Read = molden2.readlines() +v=0 +eigenv=[] +for i in range(0,len(molden2Read)): + line = molden2Read[i] + if "Ene=" in line: + eigenv.append(float(line.split()[1])) + v+=1 + if "Occup=" in line and "0.0" in line : + testValues["e+HOMO"] = eigenv[v-2] + break +molden2.close() passTest = True diff --git a/test/H2O.APMO.FCI.py b/test/H2O.APMO.FCI.py index 78ce2d6e..47af04bc 100644 --- a/test/H2O.APMO.FCI.py +++ b/test/H2O.APMO.FCI.py @@ -40,7 +40,7 @@ diffTotalEnergy = abs(refTotalEnergy - totalEnergy) diffFCIEnergy = abs(refFCIEnergy - FCIEnergy) -if (diffTotalEnergy <= 1E-8 and diffFCIEnergy <= 1E-8): +if (diffTotalEnergy <= 1E-8 and diffFCIEnergy <= 1E-6): print(testName + str_green(" ... OK")) else: print(testName + str_red(" ... NOT OK")) diff --git a/test/HCOOPs.HF.DensCube.lowdin b/test/HCOOPs.HF.DensCube.lowdin index 0626cb09..335a8d4c 100644 --- a/test/HCOOPs.HF.DensCube.lowdin +++ b/test/HCOOPs.HF.DensCube.lowdin @@ -22,7 +22,7 @@ END CONTROL OUTPUTS densityCube cubeSize=5 point1=0.0 0.0 0.0 species="E-" - densityCube cubeSize=20 point1=0.0 0.0 -2.5 species="POSITRON" + densityCube cubeSize=20 point1=0.0 0.0 -2.5 species="E+" END OUTPUTS diff --git a/test/HCOOPs.HF.DensCube.py b/test/HCOOPs.HF.DensCube.py index 43ac5465..9233a570 100644 --- a/test/HCOOPs.HF.DensCube.py +++ b/test/HCOOPs.HF.DensCube.py @@ -13,7 +13,7 @@ inputName = testName + ".lowdin" outputName = testName + ".out" cube1Name = testName + ".E-.dens.cub" -cube2Name = testName + ".POSITRON.2.dens.cub" +cube2Name = testName + ".POSITRON.dens.cub" # Reference values and tolerance refValues = { diff --git a/test/He.POTENTIALS.lowdin b/test/He.POTENTIALS.lowdin index 93ff5807..fc854ac5 100644 --- a/test/He.POTENTIALS.lowdin +++ b/test/He.POTENTIALS.lowdin @@ -36,7 +36,6 @@ END EXTERPOTENTIAL CONTROL readCoefficients=F transformToCenterOfMass=F - !CIdiagonalizationMethod="ARPACK" !CIdiagonalizationMethod="JADAMILU" CISizeOfGuessMatrix=1 convergenceMethod=1 diff --git a/test/Mg.e+.gribakin.lowdin b/test/Mg.e+.gribakin.lowdin deleted file mode 100644 index e3eba0d6..00000000 --- a/test/Mg.e+.gribakin.lowdin +++ /dev/null @@ -1,26 +0,0 @@ - -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/PsF.HF.densPlot.lowdin b/test/PsF.HF.densOrbPlots.lowdin similarity index 67% rename from test/PsF.HF.densPlot.lowdin rename to test/PsF.HF.densOrbPlots.lowdin index d9f4fd85..01c72254 100644 --- a/test/PsF.HF.densPlot.lowdin +++ b/test/PsF.HF.densOrbPlots.lowdin @@ -17,6 +17,8 @@ END CONTROL OUTPUTS densityPlot dimensions=2 point1=0.0 0.0 -10.0 point2= 0.0 0.0 10.0 + orbitalPlot dimensions=2 point1=0.0 0.0 -10.0 point2= 0.0 0.0 10.0 species="E+" orbital=1 + orbitalPlot dimensions=2 point1=0.0 0.0 -2.0 point2= 0.0 0.0 2.0 species="E-" orbital=2 END OUTPUTS diff --git a/test/PsF.HF.densOrbPlots.py b/test/PsF.HF.densOrbPlots.py new file mode 100644 index 00000000..5ee3e329 --- /dev/null +++ b/test/PsF.HF.densOrbPlots.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import sys +from colorstring import * + +if len(sys.argv)==2: + lowdinbin = sys.argv[1] +else: + lowdinbin = "lowdin2" + +testName = sys.argv[0][:-3] +inputName = testName + ".lowdin" +outputName = testName + ".out" +densplot1Name = testName + ".E-.2D.dens" +densplot2Name = testName + ".POSITRON.2D.dens" +orbplot1Name = testName + ".E-.2D.orb2" +orbplot2Name = testName + ".POSITRON.2D.orb1" +# Reference values and tolerance + +refValues = { +"HF energy" : [-99.635031860198,1E-8], +"Num e- in densplot" : [9.99830528,1E-3], +"Num e- in orbplot" : [0.99956206,1E-3], +"Num e+ in densplot" : [0.99999922,1E-5], +"Num e+ in orbplot" : [0.99999921,1E-3], +} + +testValues = dict(refValues) #copy +for value in testValues: #reset + testValues[value] = 0 #reset + +# Run calculation + +status = os.system(lowdinbin + " -i " + inputName) + +if status: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output = open(outputName, "r") +outputRead = output.readlines() + +# Values +for i in range(0,len(outputRead)): + line = outputRead[i] + if "TOTAL ENERGY =" in line: + testValues["HF energy"] = float(line.split()[3]) + +output.close() + +densplot1 = open(densplot1Name, "r") +densplot1Read = densplot1.readlines() +sumE=0 +for i in range(0,len(densplot1Read)): + line = densplot1Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[0]) + if i == 4: x2=float(values[0]) + if len(values) > 1: sumE+=float(values[0])**2*float(values[1]) +testValues["Num e- in densplot"]=2.0*3.14159265359*sumE*(x2-x1) +densplot1.close() + +densplot2 = open(densplot2Name, "r") +densplot2Read = densplot2.readlines() +sumP=0 +for i in range(0,len(densplot2Read)): + line = densplot2Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[0]) + if i == 4: x2=float(values[0]) + if len(values) > 1: sumP+=float(values[0])**2*float(values[1]) +testValues["Num e+ in densplot"]=2.0*3.14159265359*sumP*(x2-x1) +densplot2.close() + +orbplot1 = open(orbplot1Name, "r") +orbplot1Read = orbplot1.readlines() +sumE=0 +for i in range(0,len(orbplot1Read)): + line = orbplot1Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[0]) + if i == 4: x2=float(values[0]) + if len(values) > 1: sumE+=float(values[0])**2*float(values[1])**2 +testValues["Num e- in orbplot"]=2.0*3.14159265359*sumE*(x2-x1) +orbplot1.close() + +orbplot2 = open(orbplot2Name, "r") +orbplot2Read = orbplot2.readlines() +sumP=0 +for i in range(0,len(orbplot2Read)): + line = orbplot2Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[0]) + if i == 4: x2=float(values[0]) + if len(values) > 1: sumP+=float(values[0])**2*float(values[1])**2 +testValues["Num e+ in orbplot"]=2.0*3.14159265359*sumP*(x2-x1) +orbplot2.close() + +passTest = True + +for value in refValues: + diffValue = abs(refValues[value][0] - testValues[value]) + if ( diffValue <= refValues[value][1] ): + passTest = passTest * True + else : + passTest = passTest * False + print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue)) + +if passTest : + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + diff --git a/test/PsF2.CISD+molden.lowdin b/test/PsF2.CISD+molden.lowdin new file mode 100644 index 00000000..be494553 --- /dev/null +++ b/test/PsF2.CISD+molden.lowdin @@ -0,0 +1,32 @@ +GEOMETRY + e-(F) aug-cc-pvdz 0.00 0.00 -0.975 addParticles=1 multiplicity=2 + e-(F) aug-cc-pvdz 0.00 0.00 0.975 + e+ PSX-DZ 0.00 0.00 -0.975 + e+ PSX-DZ 0.00 0.00 0.975 addParticles=-1 + F dirac 0.00 0.00 -0.975 + F dirac 0.00 0.00 0.975 +END GEOMETRY + +TASKS + method = "UHF" + configurationInteractionLevel ="CISD+" +END TASKS + +CONTROL +readCoefficients=F +numberOfCIstates=1 +CINaturalOrbitals=T +CIStatesToPrint = 1 +CIdiagonalizationMethod = "JADAMILU" +totalEnergyTolerance=1E-10 +END CONTROL + +INPUT_CI + species="E-ALPHA" core=1 active=10 + species="E-BETA" core=1 active=10 + species="POSITRON" core=0 active=10 +END INPUT_CI + +OUTPUTS + moldenFile +END OUTPUTS \ No newline at end of file diff --git a/test/PsF2.CISD+molden.py b/test/PsF2.CISD+molden.py new file mode 100644 index 00000000..7a7941ba --- /dev/null +++ b/test/PsF2.CISD+molden.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import sys +from colorstring import * + +if len(sys.argv)==2: + lowdinbin = sys.argv[1] +else: + lowdinbin = "lowdin2" + +testName = sys.argv[0][:-3] +inputName = testName + ".lowdin" +outputName = testName + ".out" +molden1Name = testName + ".E-ALPHA.molden" +molden2Name = testName + ".E-BETA.molden" +molden3Name = testName + ".POSITRON.molden" +# Reference values and tolerance + +refValues = { +"HF energy" : [-198.967220107085,1E-8], +"CISD+ energy" : [-198.977666606058,1E-6], +"Natural Occ 10 e-alpha 1" : [0.9987,1E-4], +"Natural Occ 11 e-alpha 1" : [0.0017,1E-4], +"Natural Occ 9 e-beta 1" : [0.9648,1E-4], +"Natural Occ 10 e-beta 1" : [0.0344,1E-4], +"Natural Occ 1 e+ 1" : [0.9653,1E-4], +"Natural Occ 2 e+ 1" : [0.0327,1E-4], +} + +testValues = dict(refValues) #copy +for value in testValues: #reset + testValues[value] = 0 #reset + +# Run calculation + +status = os.system(lowdinbin + " -i " + inputName) + +if status: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output = open(outputName, "r") +outputRead = output.readlines() + +# Values +for i in range(0,len(outputRead)): + line = outputRead[i] + if "TOTAL ENERGY =" in line: + testValues["HF energy"] = float(line.split()[3]) + if "STATE: 1 ENERGY =" in line: + testValues["CISD+ energy"] = float(line.split()[4]) +output.close() + +molden1 = open(molden1Name, "r") +molden1Read = molden1.readlines() +v=0 +eigenv=[] +for i in range(0,len(molden1Read)): + line = molden1Read[i] + if "Occup=" in line: + v+=1 + if v==10: + testValues["Natural Occ 10 e-alpha 1"] = abs(float(line.split()[1])) + if v==11: + testValues["Natural Occ 11 e-alpha 1"] = abs(float(line.split()[1])) + if "0.0" in line: + break +molden1.close() + +molden2 = open(molden2Name, "r") +molden2Read = molden2.readlines() +v=0 +eigenv=[] +for i in range(0,len(molden2Read)): + line = molden2Read[i] + if "Occup=" in line: + v+=1 + if v==9: + testValues["Natural Occ 9 e-beta 1"] = abs(float(line.split()[1])) + if v==10: + testValues["Natural Occ 10 e-beta 1"] = abs(float(line.split()[1])) + if "0.0" in line: + break +molden2.close() + +molden3 = open(molden3Name, "r") +molden3Read = molden3.readlines() +v=0 +eigenv=[] +for i in range(0,len(molden3Read)): + line = molden3Read[i] + if "Occup=" in line: + v+=1 + if v==1: + testValues["Natural Occ 1 e+ 1"] = abs(float(line.split()[1])) + if v==2: + testValues["Natural Occ 2 e+ 1"] = abs(float(line.split()[1])) + if "0.0" in line: + break +molden3.close() + +passTest = True + +for value in refValues: + diffValue = abs(refValues[value][0] - testValues[value]) + if ( diffValue <= refValues[value][1] ): + passTest = passTest * True + else : + passTest = passTest * False + print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue)) + +if passTest : + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output.close() diff --git a/test/PsH.FCI.IT-E.lowdin b/test/PsH.FCI.IT-E.lowdin index 6b7b1e95..f8865b38 100644 --- a/test/PsH.FCI.IT-E.lowdin +++ b/test/PsH.FCI.IT-E.lowdin @@ -17,7 +17,6 @@ CINaturalOrbitals=T CIStatesToPrint = 1 !CIdiagonalizationMethod = "DSYEVX" CIdiagonalizationMethod = "JADAMILU" - !CIdiagonalizationMethod = "ARPACK" !CIPrintEigenVectorsFormat = "NONE" CIPrintEigenVectorsFormat = "OCCUPIED" !CIPrintEigenVectorsFormat = "ORBITALS" diff --git a/test/PsH.FCI.lowdin b/test/PsH.FCI.lowdin index e5b0f883..6be15f7f 100644 --- a/test/PsH.FCI.lowdin +++ b/test/PsH.FCI.lowdin @@ -17,7 +17,6 @@ CINaturalOrbitals=T CIStatesToPrint = 1 !CIdiagonalizationMethod = "DSYEVX" CIdiagonalizationMethod = "JADAMILU" - !CIdiagonalizationMethod = "ARPACK" !CIPrintEigenVectorsFormat = "NONE" CIPrintEigenVectorsFormat = "OCCUPIED" !CIPrintEigenVectorsFormat = "ORBITALS" diff --git a/test/PsOH.HF.densOrbPlots.lowdin b/test/PsOH.HF.densOrbPlots.lowdin new file mode 100644 index 00000000..0cf90678 --- /dev/null +++ b/test/PsOH.HF.densOrbPlots.lowdin @@ -0,0 +1,28 @@ + +GEOMETRY +E-(O) AUG-CC-PVTZ 0.0 0.0 0.000 addParticles=1 +E-(H) AUG-CC-PVTZ 0.0 0.0 0.971 +E+ PSX-TZ 0.0 0.0 0.000 addParticles=-1 +E+ PSX-DZ 0.0 0.0 0.971 +O dirac 0.0 0.0 0.000 +H dirac 0.0 0.0 0.971 +END GEOMETRY + +TASKS + method = "RHF" +END TASKS + +CONTROL + readCoefficients=F + numberOfPointsPerDimension=250 + totalEnergyTolerance=1E-12 +END CONTROL + +OUTPUTS + densityPlot dimensions=3 point1=0.0 -1.5 -1.5 point2= 0.0 1.5 -1.5 point3= 0.0 -1.5 1.5 species="E-" + orbitalPlot dimensions=3 point1=0.0 -2.0 -2.0 point2= 0.0 2.0 -2.0 point3= 0.0 -2.0 2.0 species="E-" orbital=3 + densityPlot dimensions=3 point1=0.0 -5.0 -5.0 point2= 0.0 5.0 -5.0 point3= 0.0 -5.0 5.0 species="E+" + orbitalPlot dimensions=3 point1=0.0 -5.0 -5.0 point2= 0.0 5.0 -5.0 point3= 0.0 -5.0 5.0 species="E+" orbital=1 +END OUTPUTS + + diff --git a/test/PsOH.HF.densOrbPlots.py b/test/PsOH.HF.densOrbPlots.py new file mode 100644 index 00000000..e79be372 --- /dev/null +++ b/test/PsOH.HF.densOrbPlots.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import sys +from colorstring import * + +if len(sys.argv)==2: + lowdinbin = sys.argv[1] +else: + lowdinbin = "lowdin2" + +testName = sys.argv[0][:-3] +inputName = testName + ".lowdin" +outputName = testName + ".out" +densplot1Name = testName + ".E-.3D.dens" +densplot2Name = testName + ".POSITRON.3D.dens" +orbplot1Name = testName + ".E-.3D.orb3" +orbplot2Name = testName + ".POSITRON.3D.orb1" +# Reference values and tolerance + +refValues = { +"HF energy" : [-75.587683637788,1E-8], +"Num e- in densplot" : [10.37978910220915,1E-1], +"Num e- in orbplot" : [1.7211264610474177,1E-1], +"Num e+ in densplot" : [0.9870409056582871,1E-2], +"Num e+ in orbplot" : [0.9870409056980444,1E-2], +} + +testValues = dict(refValues) #copy +for value in testValues: #reset + testValues[value] = 0 #reset + +# Run calculation + +status = os.system(lowdinbin + " -i " + inputName) + +if status: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output = open(outputName, "r") +outputRead = output.readlines() + +# Values +for i in range(0,len(outputRead)): + line = outputRead[i] + if "TOTAL ENERGY =" in line: + testValues["HF energy"] = float(line.split()[3]) + +output.close() + +densplot1 = open(densplot1Name, "r") +densplot1Read = densplot1.readlines() +sumE=0 +for i in range(0,len(densplot1Read)): + line = densplot1Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[1]) + if i == 4: x2=float(values[1]) + if len(values) > 1: sumE+=(float(values[0])**2+float(values[1])**2)**(1.0/2.0)*float(values[2]) +testValues["Num e- in densplot"]=2.0*sumE*(x2-x1)**2 +densplot1.close() + +densplot2 = open(densplot2Name, "r") +densplot2Read = densplot2.readlines() +sumP=0 +for i in range(0,len(densplot2Read)): + line = densplot2Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[1]) + if i == 4: x2=float(values[1]) + if len(values) > 1: sumP+=(float(values[0])**2+float(values[1])**2)**(1.0/2.0)*float(values[2]) +testValues["Num e+ in densplot"]=2.0*sumP*(x2-x1)**2 +densplot2.close() + +orbplot1 = open(orbplot1Name, "r") +orbplot1Read = orbplot1.readlines() +sumE=0 +for i in range(0,len(orbplot1Read)): + line = orbplot1Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[1]) + if i == 4: x2=float(values[1]) + if len(values) > 1: sumE+=(float(values[0])**2+float(values[1])**2)**(1.0/2.0)*float(values[2])**2 +testValues["Num e- in orbplot"]=2.0*sumE*(x2-x1)**2 +orbplot1.close() + +orbplot2 = open(orbplot2Name, "r") +orbplot2Read = orbplot2.readlines() +sumP=0 +for i in range(0,len(orbplot2Read)): + line = orbplot2Read[i] + if i > 1: + values = line.split() + if i == 3: x1=float(values[1]) + if i == 4: x2=float(values[1]) + if len(values) > 1: sumP+=(float(values[0])**2+float(values[1])**2)**(1.0/2.0)*float(values[2])**2 +testValues["Num e+ in orbplot"]=2.0*sumP*(x2-x1)**2 +orbplot2.close() + +passTest = True + +for value in refValues: + diffValue = abs(refValues[value][0] - testValues[value]) + if ( diffValue <= refValues[value][1] ): + passTest = passTest * True + else : + passTest = passTest * False + print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue)) + +if passTest : + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + diff --git a/test/PsOH.HF.densPlot.lowdin b/test/PsOH.HF.densPlot.lowdin deleted file mode 100644 index b2ed5748..00000000 --- a/test/PsOH.HF.densPlot.lowdin +++ /dev/null @@ -1,25 +0,0 @@ - -GEOMETRY -E-(O) AUG-CC-PVTZ 0.0 0.0 0.000 -E-(H) AUG-CC-PVTZ 0.0 0.0 0.971 addParticles=1 -E+ PSX-TZ 0.0 0.0 0.000 addParticles=-1 -E+ PSX-DZ 0.0 0.0 0.971 -O dirac 0.0 0.0 0.000 -H dirac 0.0 0.0 0.971 -END GEOMETRY - -TASKS - method = "RHF" -END TASKS - -CONTROL - readCoefficients=F - numberOfPointsPerDimension=250 - totalEnergyTolerance=1E-12 -END CONTROL - -OUTPUTS - densityPlot dimensions=3 point1=0.0 -10.0 -10.0 point2= 0.0 10.0 -10.0 point3= 0.0 -10.0 10.0 -END OUTPUTS - -