diff --git a/examples/basis.B.cc-pvdz b/examples/basis.B.cc-pvdz new file mode 100644 index 0000000..b9d548b --- /dev/null +++ b/examples/basis.B.cc-pvdz @@ -0,0 +1,30 @@ +1 6 +S 8 + 1 4570.0000000 0.0006960 + 2 685.9000000 0.0053530 + 3 156.5000000 0.0271340 + 4 44.4700000 0.1013800 + 5 14.4800000 0.2720550 + 6 5.1310000 0.4484030 + 7 1.8980000 0.2901230 + 8 0.3329000 0.0143220 +S 8 + 1 4570.0000000 -0.0001390 + 2 685.9000000 -0.0010970 + 3 156.5000000 -0.0054440 + 4 44.4700000 -0.0219160 + 5 14.4800000 -0.0597510 + 6 5.1310000 -0.1387320 + 7 1.8980000 -0.1314820 + 8 0.3329000 0.5395260 +S 1 + 1 0.1043000 1.0000000 +P 3 + 1 6.0010000 0.0354810 + 2 1.2410000 0.1980720 + 3 0.3364000 0.5052300 +P 1 + 1 0.0953800 1.0000000 +D 1 + 1 0.3430000 1.0000000 + diff --git a/examples/basis.Be.6-31g b/examples/basis.Be.6-31g new file mode 100644 index 0000000..dd9c9af --- /dev/null +++ b/examples/basis.Be.6-31g @@ -0,0 +1,22 @@ +1 5 +S 6 + 1 1264.5857000 0.0019448 + 2 189.9368100 0.0148351 + 3 43.1590890 0.0720906 + 4 12.0986630 0.2371542 + 5 3.8063232 0.4691987 + 6 1.2728903 0.3565202 +S 3 + 1 3.1964631 -0.1126487 + 2 0.7478133 -0.2295064 + 3 0.2199663 1.1869167 +P 3 + 1 3.1964631 0.0559802 + 2 0.7478133 0.2615506 + 3 0.2199663 0.7939723 +S 1 + 1 0.0823099 1.0000000 +P 1 + 1 0.0823099 1.0000000 + + diff --git a/examples/basis.Be.6-31g_star b/examples/basis.Be.6-31g_star new file mode 100644 index 0000000..b49a51b --- /dev/null +++ b/examples/basis.Be.6-31g_star @@ -0,0 +1,23 @@ +1 6 +S 6 + 1 1264.5857000 0.0019448 + 2 189.9368100 0.0148351 + 3 43.1590890 0.0720906 + 4 12.0986630 0.2371542 + 5 3.8063232 0.4691987 + 6 1.2728903 0.3565202 +S 3 + 1 3.1964631 -0.1126487 + 2 0.7478133 -0.2295064 + 3 0.2199663 1.1869167 +P 3 + 1 3.1964631 0.0559802 + 2 0.7478133 0.2615506 + 3 0.2199663 0.7939723 +S 1 + 1 0.0823099 1.0000000 +P 1 + 1 0.0823099 1.0000000 +D 1 + 1 0.4000000 1.0000000 + diff --git a/examples/basis.C.cc-pvdz b/examples/basis.C.cc-pvdz new file mode 100644 index 0000000..7a2f348 --- /dev/null +++ b/examples/basis.C.cc-pvdz @@ -0,0 +1,32 @@ +1 6 +S 9 +1 6.665000E+03 6.920000E-04 +2 1.000000E+03 5.329000E-03 +3 2.280000E+02 2.707700E-02 +4 6.471000E+01 1.017180E-01 +5 2.106000E+01 2.747400E-01 +6 7.495000E+00 4.485640E-01 +7 2.797000E+00 2.850740E-01 +8 5.215000E-01 1.520400E-02 +9 1.596000E-01 -3.191000E-03 +S 9 +1 6.665000E+03 -1.460000E-04 +2 1.000000E+03 -1.154000E-03 +3 2.280000E+02 -5.725000E-03 +4 6.471000E+01 -2.331200E-02 +5 2.106000E+01 -6.395500E-02 +6 7.495000E+00 -1.499810E-01 +7 2.797000E+00 -1.272620E-01 +8 5.215000E-01 5.445290E-01 +9 1.596000E-01 5.804960E-01 +S 1 +1 1.596000E-01 1.000000E+00 +P 4 +1 9.439000E+00 3.810900E-02 +2 2.002000E+00 2.094800E-01 +3 5.456000E-01 5.085570E-01 +4 1.517000E-01 4.688420E-01 +P 1 +1 1.517000E-01 1.000000E+00 +D 1 +1 5.500000E-01 1.0000000 diff --git a/examples/basis.C2.cc-pvdz b/examples/basis.C2.cc-pvdz new file mode 100644 index 0000000..16d91fa --- /dev/null +++ b/examples/basis.C2.cc-pvdz @@ -0,0 +1,64 @@ +1 6 +S 9 +1 6.665000E+03 6.920000E-04 +2 1.000000E+03 5.329000E-03 +3 2.280000E+02 2.707700E-02 +4 6.471000E+01 1.017180E-01 +5 2.106000E+01 2.747400E-01 +6 7.495000E+00 4.485640E-01 +7 2.797000E+00 2.850740E-01 +8 5.215000E-01 1.520400E-02 +9 1.596000E-01 -3.191000E-03 +S 9 +1 6.665000E+03 -1.460000E-04 +2 1.000000E+03 -1.154000E-03 +3 2.280000E+02 -5.725000E-03 +4 6.471000E+01 -2.331200E-02 +5 2.106000E+01 -6.395500E-02 +6 7.495000E+00 -1.499810E-01 +7 2.797000E+00 -1.272620E-01 +8 5.215000E-01 5.445290E-01 +9 1.596000E-01 5.804960E-01 +S 1 +1 1.596000E-01 1.000000E+00 +P 4 +1 9.439000E+00 3.810900E-02 +2 2.002000E+00 2.094800E-01 +3 5.456000E-01 5.085570E-01 +4 1.517000E-01 4.688420E-01 +P 1 +1 1.517000E-01 1.000000E+00 +D 1 +1 5.500000E-01 1.0000000 +2 6 +S 9 +1 6.665000E+03 6.920000E-04 +2 1.000000E+03 5.329000E-03 +3 2.280000E+02 2.707700E-02 +4 6.471000E+01 1.017180E-01 +5 2.106000E+01 2.747400E-01 +6 7.495000E+00 4.485640E-01 +7 2.797000E+00 2.850740E-01 +8 5.215000E-01 1.520400E-02 +9 1.596000E-01 -3.191000E-03 +S 9 +1 6.665000E+03 -1.460000E-04 +2 1.000000E+03 -1.154000E-03 +3 2.280000E+02 -5.725000E-03 +4 6.471000E+01 -2.331200E-02 +5 2.106000E+01 -6.395500E-02 +6 7.495000E+00 -1.499810E-01 +7 2.797000E+00 -1.272620E-01 +8 5.215000E-01 5.445290E-01 +9 1.596000E-01 5.804960E-01 +S 1 +1 1.596000E-01 1.000000E+00 +P 4 +1 9.439000E+00 3.810900E-02 +2 2.002000E+00 2.094800E-01 +3 5.456000E-01 5.085570E-01 +4 1.517000E-01 4.688420E-01 +P 1 +1 1.517000E-01 1.000000E+00 +D 1 +1 5.500000E-01 1.0000000 diff --git a/examples/basis.CH.aug-cc-pvdz b/examples/basis.CH.aug-cc-pvdz new file mode 100644 index 0000000..a93cea1 --- /dev/null +++ b/examples/basis.CH.aug-cc-pvdz @@ -0,0 +1,49 @@ +1 9 +S 8 + 1 6665.0000000 0.0006920 + 2 1000.0000000 0.0053290 + 3 228.0000000 0.0270770 + 4 64.7100000 0.1017180 + 5 21.0600000 0.2747400 + 6 7.4950000 0.4485640 + 7 2.7970000 0.2850740 + 8 0.5215000 0.0152040 +S 8 + 1 6665.0000000 -0.0001460 + 2 1000.0000000 -0.0011540 + 3 228.0000000 -0.0057250 + 4 64.7100000 -0.0233120 + 5 21.0600000 -0.0639550 + 6 7.4950000 -0.1499810 + 7 2.7970000 -0.1272620 + 8 0.5215000 0.5445290 +S 1 + 1 0.1596000 1.0000000 +S 1 + 1 0.0469000 1.0000000 +P 3 + 1 9.4390000 0.0381090 + 2 2.0020000 0.2094800 + 3 0.5456000 0.5085570 +P 1 + 1 0.1517000 1.0000000 +P 1 + 1 0.0404100 1.0000000 +D 1 + 1 0.5500000 1.0000000 +D 1 + 1 0.1510000 1.0000000 +2 5 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 +S 1 + 1 0.1220000 1.0000000 +S 1 + 1 0.0297400 1.0000000 +P 1 + 1 0.7270000 1.0000000 +P 1 + 1 0.1410000 1.0000000 + diff --git a/examples/basis.CH.cc-pvdz b/examples/basis.CH.cc-pvdz new file mode 100644 index 0000000..c7ec793 --- /dev/null +++ b/examples/basis.CH.cc-pvdz @@ -0,0 +1,38 @@ +1 6 +S 8 + 1 6665.0000000 0.0006920 + 2 1000.0000000 0.0053290 + 3 228.0000000 0.0270770 + 4 64.7100000 0.1017180 + 5 21.0600000 0.2747400 + 6 7.4950000 0.4485640 + 7 2.7970000 0.2850740 + 8 0.5215000 0.0152040 +S 8 + 1 6665.0000000 -0.0001460 + 2 1000.0000000 -0.0011540 + 3 228.0000000 -0.0057250 + 4 64.7100000 -0.0233120 + 5 21.0600000 -0.0639550 + 6 7.4950000 -0.1499810 + 7 2.7970000 -0.1272620 + 8 0.5215000 0.5445290 +S 1 + 1 0.1596000 1.0000000 +P 3 + 1 9.4390000 0.0381090 + 2 2.0020000 0.2094800 + 3 0.5456000 0.5085570 +P 1 + 1 0.1517000 1.0000000 +D 1 + 1 0.5500000 1.0000000 +2 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 +S 1 + 1 0.1220000 1.0000000 +P 1 + 1 0.7270000 1.0000000 diff --git a/examples/basis.CO.sto-3g b/examples/basis.CO.sto-3g new file mode 100644 index 0000000..d69bfd4 --- /dev/null +++ b/examples/basis.CO.sto-3g @@ -0,0 +1,28 @@ +1 3 +S 3 + 1 71.6168370 0.15432897 + 2 13.0450960 0.53532814 + 3 3.5305122 0.44463454 +S 3 + 1 2.9412494 -0.0999672 + 2 0.6834831 0.3995128 + 3 0.2222899 0.7001155 +P 3 + 1 2.9412494 0.1559163 + 2 0.6834831 0.6076837 + 3 0.2222899 0.3919574 +2 3 +S 3 + 1 130.7093200 0.15432897 + 2 23.8088610 0.53532814 + 3 6.4436083 0.44463454 +S 3 + 1 5.0331513 -0.0999672 + 2 1.1695961 0.3995128 + 3 0.3803890 0.7001155 +P 3 + 1 5.0331513 0.1559163 + 2 1.1695961 0.6076837 + 3 0.3803890 0.3919574 + + diff --git a/examples/basis.F.cc-pvdz b/examples/basis.F.cc-pvdz new file mode 100644 index 0000000..82e0afb --- /dev/null +++ b/examples/basis.F.cc-pvdz @@ -0,0 +1,29 @@ +1 6 +S 8 + 1 14710.0000000 0.0007210 + 2 2207.0000000 0.0055530 + 3 502.8000000 0.0282670 + 4 142.6000000 0.1064440 + 5 46.4700000 0.2868140 + 6 16.7000000 0.4486410 + 7 6.3560000 0.2647610 + 8 1.3160000 0.0153330 +S 8 + 1 14710.0000000 -0.0001650 + 2 2207.0000000 -0.0013080 + 3 502.8000000 -0.0064950 + 4 142.6000000 -0.0266910 + 5 46.4700000 -0.0736900 + 6 16.7000000 -0.1707760 + 7 6.3560000 -0.1123270 + 8 1.3160000 0.5628140 +S 1 + 1 0.3897000 1.0000000 +P 3 + 1 22.6700000 0.0448780 + 2 4.9770000 0.2357180 + 3 1.3470000 0.5085210 +P 1 + 1 0.3471000 1.0000000 +D 1 + 1 1.6400000 1.0000000 diff --git a/examples/basis.H.cc-pvdz b/examples/basis.H.cc-pvdz new file mode 100644 index 0000000..79a747a --- /dev/null +++ b/examples/basis.H.cc-pvdz @@ -0,0 +1,9 @@ +1 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 +S 1 + 1 0.1220000 1.0000000 +P 1 + 1 0.7270000 1.0000000 diff --git a/examples/basis.H.sto-3g b/examples/basis.H.sto-3g new file mode 100644 index 0000000..10ed3c2 --- /dev/null +++ b/examples/basis.H.sto-3g @@ -0,0 +1,5 @@ +1 1 +S 3 + 1 3.42525091 0.15432897 + 2 0.62391373 0.53532814 + 3 0.16885540 0.44463454 diff --git a/examples/basis.H2.6-311g b/examples/basis.H2.6-311g new file mode 100644 index 0000000..2ca2f7c --- /dev/null +++ b/examples/basis.H2.6-311g @@ -0,0 +1,18 @@ +1 3 +S 3 + 1 33.8650000 0.0254938 + 2 5.0947900 0.1903730 + 3 1.1587900 0.8521610 +S 1 + 1 0.3258400 1.0000000 +S 1 + 1 0.1027410 1.0000000 +2 3 +S 3 + 1 33.8650000 0.0254938 + 2 5.0947900 0.1903730 + 3 1.1587900 0.8521610 +S 1 + 1 0.3258400 1.0000000 +S 1 + 1 0.1027410 1.0000000 diff --git a/examples/basis.H2.pipis b/examples/basis.H2.pipis new file mode 100644 index 0000000..d125feb --- /dev/null +++ b/examples/basis.H2.pipis @@ -0,0 +1,6 @@ +1 1 +P 1 + 1 0.7270000 1.0000000 +2 1 +P 1 + 1 0.7270000 1.0000000 diff --git a/examples/basis.HeH+.sto-3g b/examples/basis.HeH+.sto-3g new file mode 100644 index 0000000..a21d2b0 --- /dev/null +++ b/examples/basis.HeH+.sto-3g @@ -0,0 +1,11 @@ +1 1 +S 3 + 1 6.36242139 0.15432897 + 2 1.15892300 0.53532814 + 3 0.31364979 0.44463454 +2 1 +S 3 + 1 3.42525091 0.15432897 + 2 0.62391373 0.53532814 + 3 0.16885540 0.44463454 + diff --git a/examples/basis.Li+.cc-pvdz b/examples/basis.Li+.cc-pvdz new file mode 100644 index 0000000..b2b2293 --- /dev/null +++ b/examples/basis.Li+.cc-pvdz @@ -0,0 +1,30 @@ +1 6 +S 8 + 1 1469.0000000 0.0007660 + 2 220.5000000 0.0058920 + 3 50.2600000 0.0296710 + 4 14.2400000 0.1091800 + 5 4.5810000 0.2827890 + 6 1.5800000 0.4531230 + 7 0.5640000 0.2747740 + 8 0.0734500 0.0097510 +S 8 + 1 1469.0000000 -0.0001200 + 2 220.5000000 -0.0009230 + 3 50.2600000 -0.0046890 + 4 14.2400000 -0.0176820 + 5 4.5810000 -0.0489020 + 6 1.5800000 -0.0960090 + 7 0.5640000 -0.1363800 + 8 0.0734500 0.5751020 +S 1 + 1 0.0280500 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 +P 1 + 1 0.0240300 1.0000000 +D 1 + 1 0.1239000 1.0000000 + diff --git a/examples/basis.Li.cc-pvdz b/examples/basis.Li.cc-pvdz new file mode 100644 index 0000000..b2b2293 --- /dev/null +++ b/examples/basis.Li.cc-pvdz @@ -0,0 +1,30 @@ +1 6 +S 8 + 1 1469.0000000 0.0007660 + 2 220.5000000 0.0058920 + 3 50.2600000 0.0296710 + 4 14.2400000 0.1091800 + 5 4.5810000 0.2827890 + 6 1.5800000 0.4531230 + 7 0.5640000 0.2747740 + 8 0.0734500 0.0097510 +S 8 + 1 1469.0000000 -0.0001200 + 2 220.5000000 -0.0009230 + 3 50.2600000 -0.0046890 + 4 14.2400000 -0.0176820 + 5 4.5810000 -0.0489020 + 6 1.5800000 -0.0960090 + 7 0.5640000 -0.1363800 + 8 0.0734500 0.5751020 +S 1 + 1 0.0280500 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 +P 1 + 1 0.0240300 1.0000000 +D 1 + 1 0.1239000 1.0000000 + diff --git a/examples/basis.Li.sto-3g b/examples/basis.Li.sto-3g new file mode 100644 index 0000000..f03c9ff --- /dev/null +++ b/examples/basis.Li.sto-3g @@ -0,0 +1,14 @@ +1 3 +S 3 + 1 16.1195750 0.15432897 + 2 2.9362007 0.53532814 + 3 0.7946505 0.44463454 +S 3 + 1 0.6362897 -0.0999672 + 2 0.1478601 0.3995128 + 3 0.0480887 0.7001155 +P 3 + 1 0.6362897 0.1559163 + 2 0.1478601 0.6076837 + 3 0.0480887 0.3919574 + diff --git a/examples/basis.Ne.6-31g b/examples/basis.Ne.6-31g new file mode 100644 index 0000000..f1a1e79 --- /dev/null +++ b/examples/basis.Ne.6-31g @@ -0,0 +1,21 @@ +1 5 +S 6 + 1 8425.8515300 0.0018843481 + 2 1268.5194000 0.0143368994 + 3 289.6214140 0.0701096233 + 4 81.8590040 0.2373732660 + 5 26.2515079 0.4730071260 + 6 9.09472051 0.3484012410 +S 3 + 1 26.5321310 -0.1071183 + 2 6.1017550 -0.1461638 + 3 1.6962715 1.1277735 +P 3 + 1 26.5321310 0.0719096 + 2 6.1017550 0.3495134 + 3 1.6962715 0.7199405 +S 1 + 1 0.4458187 1.0000000 +P 1 + 1 0.4458187 1.0000000 + diff --git a/examples/basis.Ne.6-31g_star b/examples/basis.Ne.6-31g_star new file mode 100644 index 0000000..4fe8cec --- /dev/null +++ b/examples/basis.Ne.6-31g_star @@ -0,0 +1,23 @@ +1 6 +S 6 + 1 8425.8515300 0.0018843481 + 2 1268.5194000 0.0143368994 + 3 289.6214140 0.0701096233 + 4 81.8590040 0.2373732660 + 5 26.2515079 0.4730071260 + 6 9.09472051 0.3484012410 +S 3 + 1 26.5321310 -0.1071183 + 2 6.1017550 -0.1461638 + 3 1.6962715 1.1277735 +P 3 + 1 26.5321310 0.0719096 + 2 6.1017550 0.3495134 + 3 1.6962715 0.7199405 +S 1 + 1 0.4458187 1.0000000 +P 1 + 1 0.4458187 1.0000000 +D 1 + 1 0.8000000 1.0000000 + diff --git a/examples/molecule.B b/examples/molecule.B new file mode 100644 index 0000000..77f6e05 --- /dev/null +++ b/examples/molecule.B @@ -0,0 +1,4 @@ +# nAt nEla nElb nCore nRyd + 1 3 2 0 0 +# Znuc x y z + B 0. 0. 0. diff --git a/examples/molecule.C b/examples/molecule.C new file mode 100644 index 0000000..7f71f46 --- /dev/null +++ b/examples/molecule.C @@ -0,0 +1,4 @@ +# nAt nEla nElb nCore nRyd + 1 5 1 0 0 +# Znuc x y z + C 0. 0. 0. diff --git a/examples/molecule.CH b/examples/molecule.CH new file mode 100644 index 0000000..e2a4fd3 --- /dev/null +++ b/examples/molecule.CH @@ -0,0 +1,5 @@ +# nAt nEla nElb nCore nRyd + 2 4 3 0 0 +# Znuc x y z + C 0. 0. -0.16245872 + H 0. 0. 1.93436816 diff --git a/examples/molecule.F b/examples/molecule.F new file mode 100644 index 0000000..52f335d --- /dev/null +++ b/examples/molecule.F @@ -0,0 +1,4 @@ +# nAt nEla nElb nCore nRyd + 1 5 4 0 0 +# Znuc x y z + F 0. 0. 0. diff --git a/examples/molecule.H b/examples/molecule.H new file mode 100644 index 0000000..fd4bfbe --- /dev/null +++ b/examples/molecule.H @@ -0,0 +1,4 @@ +# nAt nEla nElb nCore nRyd + 1 1 0 0 0 +# Znuc x y z + H 0. 0. 0. diff --git a/examples/molecule.Li+ b/examples/molecule.Li+ new file mode 100644 index 0000000..d8493e9 --- /dev/null +++ b/examples/molecule.Li+ @@ -0,0 +1,4 @@ +# nAt nEla nElb nCore nRyd + 1 1 1 0 0 +# Znuc x y z + Li 0.0 0.0 0.0 diff --git a/input/methods b/input/methods index 61f2d19..f7bba22 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - T F F + F T F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) @@ -7,9 +7,9 @@ # drCCD rCCD lCCD pCCD F F F F # CIS* CIS(D) CID CISD - F F F F + T F F F # RPA* RPAx* ppRPA - F T F + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW diff --git a/input/options b/input/options index 4b6f510..4dfc2da 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T F T + T T T T F # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index cc06f85..b846005 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -1,5 +1,5 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & - dipole_int_aa,dipole_int_bb,eHF) + dipole_int_aa,dipole_int_bb,eHF,cHF,S) ! Perform configuration interaction single calculation` @@ -17,6 +17,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) @@ -27,7 +29,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E ! Local variables logical :: dump_matrix = .false. - logical :: dump_trans = .false. + logical :: dump_trans = .false. integer :: ispin double precision :: lambda @@ -80,8 +82,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E call diagonalize_matrix(nS_sc,A_sc,Omega_sc) call print_excitation('UCIS ',5,nS_sc,Omega_sc) - call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - Omega_sc,transpose(A_sc),transpose(A_sc)) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + cHF,S,Omega_sc,transpose(A_sc),transpose(A_sc)) if(dump_trans) then print*,'Spin-conserved CIS transition vectors' @@ -118,14 +120,14 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E write(*,*) endif -! call diagonalize_matrix(nS_sf,A_sf,Omega_sf) - allocate(order(nS_sf)) - call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf) - call quick_sort(Omega_sf,order(:),nS_sf) + call diagonalize_matrix(nS_sf,A_sf,Omega_sf) +! allocate(order(nS_sf)) +! call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf) +! call quick_sort(Omega_sf,order(:),nS_sf) call print_excitation('UCIS ',6,nS_sf,Omega_sf) - call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & - Omega_sf,transpose(A_sf),transpose(A_sf)) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & + cHF,S,Omega_sf,transpose(A_sf),transpose(A_sf)) if(dump_trans) then print*,'Spin-flip CIS transition vectors' diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 064d13c..67a5d8e 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -1,4 +1,4 @@ -subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) +subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) ! Print one- and two-electron energies and other stuff for UHF calculation @@ -7,7 +7,7 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) integer,intent(in) :: nBas integer,intent(in) :: nO(nspin) - double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: Ov(nBas,nBas) double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: c(nBas,nBas,nspin) double precision,intent(in) :: ENuc @@ -19,14 +19,12 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) double precision,intent(in) :: dipole(ncart) integer :: ixyz - integer :: i,j integer :: ispin double precision :: HOMO(nspin) double precision :: LUMO(nspin) double precision :: Gap(nspin) - double precision :: S2_exact - double precision :: S2 - integer :: spin_state + double precision :: S_exact,S2_exact + double precision :: S,S2 ! HOMO and LUMO @@ -47,9 +45,10 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) end do S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) - S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(S,c(:,1:nO(2),2)))**2) + S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2) - spin_state = nO(1) - nO(2) + 1 + S_exact = 0.5d0*dble(nO(1) - nO(2)) + S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2) ! Dump results @@ -91,7 +90,8 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,I6)') ' 2S+1 :',spin_state + write(*,'(A40,F13.6)') ' S (exact) :',2d0*S_exact + 1d0 + write(*,'(A40,F13.6)') ' S :',2d0*S + 1d0 write(*,'(A40,F13.6)') ' (exact) :',S2_exact write(*,'(A40,F13.6)') ' :',S2 write(*,'(A60)') '-------------------------------------------------' diff --git a/src/LR/linear_response.f90 b/src/LR/linear_response.f90 index e443d52..427b4b8 100644 --- a/src/LR/linear_response.f90 +++ b/src/LR/linear_response.f90 @@ -69,9 +69,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - do ia=1,nS - if(Omega(ia) < 0d0) Omega(ia) = 0d0 - end do +! do ia=1,nS +! if(Omega(ia) < 0d0) Omega(ia) = 0d0 +! end do call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq) call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv) @@ -83,9 +83,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') - do ia=1,nS - if(Omega(ia) < 0d0) Omega(ia) = 0d0 - end do +! do ia=1,nS +! if(Omega(ia) < 0d0) Omega(ia) = 0d0 +! end do Omega = sqrt(Omega) diff --git a/src/LR/oscillator_strength.f90 b/src/LR/oscillator_strength.f90 index 5524169..8481d1f 100644 --- a/src/LR/oscillator_strength.f90 +++ b/src/LR/oscillator_strength.f90 @@ -1,4 +1,4 @@ -subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) +subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os) ! Compute linear response @@ -13,6 +13,7 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS + integer,intent(in) :: maxS double precision :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: XpY(nS,nS) @@ -20,7 +21,6 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) ! Local variables - logical :: debug = .false. integer :: ia,jb,i,j,a,b integer :: ixyz @@ -32,7 +32,7 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) ! Memory allocation - allocate(f(nS,ncart)) + allocate(f(maxS,ncart)) ! Initialization @@ -40,7 +40,7 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) ! Compute dipole moments and oscillator strengths - do ia=1,nS + do ia=1,maxS do ixyz=1,ncart jb = 0 do j=nC+1,nO @@ -53,24 +53,19 @@ subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) end do f(:,:) = sqrt(2d0)*f(:,:) - do ia=1,nS + do ia=1,maxS os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) end do - - if(debug) then - - write(*,*) '------------------------' - write(*,*) ' Dipole moments (X Y Z) ' - write(*,*) '------------------------' - call matout(nS,ncart,f) - write(*,*) - - write(*,*) '----------------------' - write(*,*) ' Oscillator strengths ' - write(*,*) '----------------------' - call matout(nS,1,os) - write(*,*) - - end if + + write(*,*) '---------------------------------------------------------------' + write(*,*) ' Transition dipole moment (au) ' + write(*,*) '---------------------------------------------------------------' + write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.' + write(*,*) '---------------------------------------------------------------' + do ia=1,maxS + write(*,'(I3,5F12.6)') ia,(f(ia,ixyz),ixyz=1,ncart),sum(f(ia,:)**2),os(ia) + end do + write(*,*) '---------------------------------------------------------------' + write(*,*) end subroutine oscillator_strength diff --git a/src/LR/print_transition_vectors.f90 b/src/LR/print_transition_vectors.f90 index c07e4a0..c861b96 100644 --- a/src/LR/print_transition_vectors.f90 +++ b/src/LR/print_transition_vectors.f90 @@ -21,10 +21,8 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, ! Local variables - logical :: debug = .false. - integer :: ia,jb,i,j,a,b - integer :: ixyz - integer,parameter :: maxS = 10 + integer :: ia,jb,j,b + integer :: maxS = 10 double precision :: S2 double precision,parameter :: thres_vec = 0.1d0 double precision,allocatable :: X(:) @@ -33,20 +31,29 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, ! Memory allocation - allocate(X(nS),Y(nS),os(nS)) + maxS = min(nS,maxS) + allocate(X(nS),Y(nS),os(maxS)) ! Compute oscillator strengths os(:) = 0d0 - if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) + if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os) ! Print details about excitations - do ia=1,min(nS,maxS) + do ia=1,maxS X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + ! values + + if(spin_allowed) then + S2 = 0d0 + else + S2 = 2d0 + end if + print*,'-------------------------------------------------------------' write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2 diff --git a/src/LR/print_unrestricted_transition_vectors.f90 b/src/LR/print_unrestricted_transition_vectors.f90 index 11b5fb5..515248c 100644 --- a/src/LR/print_unrestricted_transition_vectors.f90 +++ b/src/LR/print_unrestricted_transition_vectors.f90 @@ -1,5 +1,5 @@ -subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, & - Omega,XpY,XmY) +subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, & + c,S,Omega,XpY,XmY) ! Print transition vectors for linear response calculation @@ -8,7 +8,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n ! Input variables - logical,intent(in) :: spin_allowed + integer,intent(in) :: ispin integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -20,94 +20,148 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n integer,intent(in) :: nSt double precision :: dipole_int_aa(nBas,nBas,ncart) double precision :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: Omega(nSt) double precision,intent(in) :: XpY(nSt,nSt) double precision,intent(in) :: XmY(nSt,nSt) ! Local variables - logical :: debug = .false. - integer :: ia,jb,i,j,a,b - integer :: ixyz - integer :: ispin - integer,parameter :: maxS = 10 - double precision :: S2 + integer :: ia,jb,j,b + integer :: maxS = 10 double precision,parameter :: thres_vec = 0.1d0 double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: f(:,:) double precision,allocatable :: os(:) + double precision,allocatable :: S2(:) ! Memory allocation - allocate(X(nSt),Y(nSt),os(nSt)) + maxS = min(nSt,maxS) + allocate(X(nSt),Y(nSt),os(maxS),S2(maxS)) ! Compute oscillator strengths os(:) = 0d0 - if(spin_allowed) call unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, & + if(ispin == 1) call unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS, & dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os) -! Print details about excitations +! Compute - do ia=1,min(nSt,maxS) + call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2) - X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) - Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) +! Print details about spin-conserved excitations - S2 = (nO(1) - nO(2))/2d0 - S2 = 2d0*S2+1d0 - S2 = 0.0d0 - do jb=1,nSa - S2 = S2 + 4d0*(X(jb)**2 + Y(jb)**2) - end do - do jb=1,nSb - S2 = S2 - 4d0*(X(nSa+jb)**2 + Y(nSa+jb)**2) - end do + if(ispin == 1) then - print*,'-------------------------------------------------------------' - write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2 - print*,'-------------------------------------------------------------' + do ia=1,maxS - ! Spin-up transitions + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb) + + print*,'-------------------------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & + ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2(ia) + print*,'-------------------------------------------------------------' + + ! Spin-up transitions + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb) + end do end do - end do - - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb) + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb) + end do end do + + ! Spin-down transitions + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(nSa+jb) + end do + end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(nSa+jb) + end do + end do + write(*,*) + end do - ! Spin-down transitions + end if - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(nSa+jb) - end do - end do - - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(nSa+jb) - end do - end do - write(*,*) +! Print details about spin-flip excitations - end do + if(ispin == 2) then + + do ia=1,maxS + + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + + print*,'-------------------------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & + ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2(ia) + print*,'-------------------------------------------------------------' + + ! Spin-up transitions + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'B = ',X(jb) + end do + end do + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'B = ',Y(jb) + end do + end do + + ! Spin-down transitions + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'B = ',X(nSa+jb) + end do + end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'B = ',Y(nSa+jb) + end do + end do + write(*,*) + + end do + + end if ! Thomas-Reiche-Kuhn sum rule diff --git a/src/LR/unrestricted_S2_expval.f90 b/src/LR/unrestricted_S2_expval.f90 new file mode 100644 index 0000000..559d6fc --- /dev/null +++ b/src/LR/unrestricted_S2_expval.f90 @@ -0,0 +1,177 @@ +subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2) + +! Compute for linear response excited states + + implicit none + include 'parameters.h' + +! Input variables + + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + integer,intent(in) :: maxS + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: XpY(nSt,nSt) + double precision,intent(in) :: XmY(nSt,nSt) + +! Local variables + + integer :: m + integer :: ia,i,a + double precision :: S2_exact + double precision :: S2_gs + double precision,allocatable :: Xa(:,:), Xb(:,:), Ya(:,:), Yb(:,:) + double precision,allocatable :: Xat(:,:),Xbt(:,:),Yat(:,:),Ybt(:,:) + double precision,allocatable :: OO(:,:), OV(:,:), VO(:,:), VV(:,:) + double precision,allocatable :: OOt(:,:),OVt(:,:),VOt(:,:),VVt(:,:) + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: S2(maxS) + +! Memory allocation + + allocate(OO(nO(1)-nC(1),nO(2)-nC(2)), OV(nO(1)-nC(1),nV(2)-nR(2)), VO(nV(1)-nR(1),nO(2)-nC(2)), VV(nV(1)-nR(1),nV(2)-nR(2)), & + OOt(nO(2)-nC(2),nO(1)-nC(1)),OVt(nV(2)-nR(2),nO(1)-nC(1)),VOt(nO(2)-nC(2),nV(1)-nR(1)),VVt(nV(2)-nR(2),nV(1)-nR(1))) + +! Overlap matrix between spin-up and spin-down orbitals + + OO(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nC(2)+1:nO(2) ,2))) + OV(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2))) + VO(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nC(2)+1:nO(2) ,2))) + VV(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2))) + + OOt(:,:) = transpose(OO(:,:)) + OVt(:,:) = transpose(OV(:,:)) + VOt(:,:) = transpose(VO(:,:)) + VVt(:,:) = transpose(VV(:,:)) + +!-------------------------! +! for ground state ! +!-------------------------! + + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2_gs = S2_exact + nO(2) - sum(OO(:,:)**2) + +!------------------------------------------! +! for spin-conserved-excited states ! +!------------------------------------------! + + if(ispin == 1) then + + allocate(Xa(nO(1)-nC(1),nV(1)-nR(1)), Ya(nO(1)-nC(1),nV(1)-nR(1)), Xb(nO(2)-nC(2),nV(2)-nR(2)), Yb(nO(2)-nC(2),nV(2)-nR(2)), & + Xat(nV(1)-nR(1),nO(1)-nC(1)),Yat(nV(1)-nR(1),nO(1)-nC(1)),Xbt(nV(2)-nR(2),nO(2)-nC(2)),Ybt(nV(2)-nR(2),nO(2)-nC(2))) + + do m=1,maxS + + ia = 0 + do i=nC(1)+1,nO(1) + do a=1,nV(1)-nR(1) + ia = ia + 1 + Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia)) + Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia)) + end do + end do + + ia = 0 + do i=nC(2)+1,nO(2) + do a=1,nV(2)-nR(2) + ia = ia + 1 + Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia)) + Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia)) + end do + end do + + Xat(:,:) = transpose(Xa(:,:)) + Xbt(:,:) = transpose(Xb(:,:)) + Yat(:,:) = transpose(Ya(:,:)) + Ybt(:,:) = transpose(Yb(:,:)) + + S2(m) = S2_gs & + + trace_matrix(nV(1),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) & + + trace_matrix(nV(2),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) & + - trace_matrix(nO(1),matmul(Xa,matmul(VO,matmul(VOt,Xat)))) & + - trace_matrix(nO(2),matmul(Xb,matmul(OVt,matmul(OV,Xbt)))) & + - 2d0*trace_matrix(nO(1),matmul(OO,matmul(Xb,matmul(VVt,Xat)))) & + + - 2d0*trace_matrix(nV(2),matmul(OVt,matmul(Xa,matmul(VO,Yb)))) & + - 2d0*trace_matrix(nV(1),matmul(VO,matmul(Xb,matmul(OVt,Ya)))) & + + - trace_matrix(nV(1),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) & + - trace_matrix(nV(2),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) & + + trace_matrix(nO(1),matmul(Ya,matmul(VO,matmul(VOt,Yat)))) & + + trace_matrix(nO(2),matmul(Yb,matmul(OVt,matmul(OV,Ybt)))) & + + 2d0*trace_matrix(nO(1),matmul(Ya,matmul(VV,matmul(Ybt,OOt)))) + + end do + + end if + +!------------------------------------------! +! for spin-conserved-excited states ! +!------------------------------------------! + + if(ispin == 2) then + + allocate(Xa(nO(1)-nC(1),nV(2)-nR(2)), Ya(nO(1)-nC(1),nV(2)-nR(2)), Xb(nO(2)-nC(2),nV(1)-nR(1)), Yb(nO(2)-nC(2),nV(1)-nR(1)), & + Xat(nV(2)-nR(2),nO(1)-nC(1)),Yat(nV(2)-nR(2),nO(1)-nC(1)),Xbt(nV(1)-nR(1),nO(2)-nC(2)),Ybt(nV(1)-nR(1),nO(2)-nC(2))) + + do m=1,maxS + + ia = 0 + do i=nC(1)+1,nO(1) + do a=1,nV(2)-nR(2) + ia = ia + 1 + Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia)) + Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia)) + end do + end do + + ia = 0 + do i=nC(2)+1,nO(2) + do a=1,nV(1)-nR(1) + ia = ia + 1 + Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia)) + Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia)) + end do + end do + + Xat(:,:) = transpose(Xa(:,:)) + Xbt(:,:) = transpose(Xb(:,:)) + Yat(:,:) = transpose(Ya(:,:)) + Ybt(:,:) = transpose(Yb(:,:)) + + S2(m) = 0d0 + +! S2(m) = S2_gs & +! + trace_matrix(nV(1),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) & +! + trace_matrix(nV(2),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) & +! - trace_matrix(nO(1),matmul(Xa,matmul(VO,matmul(VOt,Xat)))) & +! - trace_matrix(nO(2),matmul(Xb,matmul(OVt,matmul(OV,Xbt)))) & +! - 2d0*trace_matrix(nO(1),matmul(OO,matmul(Xb,matmul(VVt,Xat)))) & +! +! - 2d0*trace_matrix(nV(2),matmul(OVt,matmul(Xa,matmul(VO,Yb)))) & +! - 2d0*trace_matrix(nV(1),matmul(VO,matmul(Xb,matmul(OVt,Ya)))) & +! +! - trace_matrix(nV(1),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) & +! - trace_matrix(nV(2),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) & +! + trace_matrix(nO(1),matmul(Ya,matmul(VO,matmul(VOt,Yat)))) & +! + trace_matrix(nO(2),matmul(Yb,matmul(OVt,matmul(OV,Ybt)))) & +! + 2d0*trace_matrix(nO(1),matmul(Ya,matmul(VV,matmul(Ybt,OOt)))) + + end do + + end if + +end subroutine unrestricted_S2_expval diff --git a/src/LR/unrestricted_linear_response.f90 b/src/LR/unrestricted_linear_response.f90 index 422a863..3c747f8 100644 --- a/src/LR/unrestricted_linear_response.f90 +++ b/src/LR/unrestricted_linear_response.f90 @@ -35,7 +35,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Local variables integer :: ia - double precision :: trace_matrix + double precision,external :: trace_matrix double precision,allocatable :: A(:,:) double precision,allocatable :: B(:,:) double precision,allocatable :: ApB(:,:) @@ -90,9 +90,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - do ia=1,nSt - if(Omega(ia) < 0d0) Omega(ia) = 0d0 - end do +! do ia=1,nSt +! if(Omega(ia) < 0d0) Omega(ia) = 0d0 +! end do call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq) call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv) @@ -104,9 +104,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') - do ia=1,nSt - if(Omega(ia) < 0d0) Omega(ia) = 0d0 - end do +! do ia=1,nSt +! if(Omega(ia) < 0d0) Omega(ia) = 0d0 +! end do Omega = sqrt(Omega) diff --git a/src/LR/unrestricted_linear_response_A_matrix.f90 b/src/LR/unrestricted_linear_response_A_matrix.f90 index 6a49cc9..7878f9c 100644 --- a/src/LR/unrestricted_linear_response_A_matrix.f90 +++ b/src/LR/unrestricted_linear_response_A_matrix.f90 @@ -121,9 +121,6 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do - print*,'nSa,nSb,nSt',nSa,nSb,nSt - call matout(nSt,nSt,A_lr) - end if !----------------------------------------------- @@ -144,10 +141,8 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do j=nC(1)+1,nO(1) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 -! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B)' A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - - (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i) - + - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,b,j,a) end do end do end do @@ -165,17 +160,13 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa jb = jb + 1 A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) -! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb) + - (1d0 - delta_dRPA)*lambda*ERI_aabb(b,i,a,j) end do end do end do end do - print*,'nSa,nSb,nSt',nSa,nSb,nSt - call matout(nSt,nSt,A_lr) - end if diff --git a/src/LR/unrestricted_linear_response_B_matrix.f90 b/src/LR/unrestricted_linear_response_B_matrix.f90 index 5cf001c..d516743 100644 --- a/src/LR/unrestricted_linear_response_B_matrix.f90 +++ b/src/LR/unrestricted_linear_response_B_matrix.f90 @@ -139,7 +139,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(i,a,b,j) + B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(i,j,b,a) end do end do @@ -157,7 +157,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a) + B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_aabb(j,i,a,b) end do end do diff --git a/src/LR/unrestricted_oscillator_strength.f90 b/src/LR/unrestricted_oscillator_strength.f90 index 4a4f8b4..50d18c1 100644 --- a/src/LR/unrestricted_oscillator_strength.f90 +++ b/src/LR/unrestricted_oscillator_strength.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os) +subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os) ! Compute linear response @@ -17,6 +17,7 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt + integer,intent(in) :: maxS double precision :: dipole_int_aa(nBas,nBas,ncart) double precision :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: Omega(nSt) @@ -25,7 +26,6 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo ! Local variables - logical :: debug = .false. integer :: ia,jb,i,j,a,b integer :: ixyz @@ -33,11 +33,11 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo ! Output variables - double precision :: os(nSt) + double precision :: os(maxS) ! Memory allocation - allocate(f(nSt,ncart)) + allocate(f(maxS,ncart)) ! Initialization @@ -45,7 +45,7 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo ! Compute dipole moments and oscillator strengths - do ia=1,nSt + do ia=1,maxS do ixyz=1,ncart jb = 0 @@ -67,24 +67,19 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipo end do end do - do ia=1,nSt + do ia=1,maxS os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) end do - if(debug) then - - write(*,*) '------------------------' - write(*,*) ' Dipole moments (X Y Z) ' - write(*,*) '------------------------' - call matout(nS,ncart,f) - write(*,*) - - write(*,*) '----------------------' - write(*,*) ' Oscillator strengths ' - write(*,*) '----------------------' - call matout(nS,1,os) - write(*,*) - - end if + write(*,*) '---------------------------------------------------------------' + write(*,*) ' Transition dipole moment (au) ' + write(*,*) '---------------------------------------------------------------' + write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.' + write(*,*) '---------------------------------------------------------------' + do ia=1,maxS + write(*,'(I3,5F12.6)') ia,(f(ia,ixyz),ixyz=1,ncart),sum(f(ia,:)**2),os(ia) + end do + write(*,*) '---------------------------------------------------------------' + write(*,*) end subroutine unrestricted_oscillator_strength diff --git a/src/MBPT/unrestricted_Bethe_Salpeter.f90 b/src/MBPT/unrestricted_Bethe_Salpeter.f90 index 7e4dc02..d61182e 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter.f90 @@ -98,8 +98,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), & OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) - call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) +! call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & +! OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -131,6 +131,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) + ! Compute spin-flip BSE excitation energies + call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), & OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 index 25e81c5..31a6ef3 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -117,7 +117,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps enddo - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 2d0*lambda*chi + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi end do end do @@ -141,7 +141,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps enddo - A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 2d0*lambda*chi + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi end do end do diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 index 2f559ee..69965f2 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 @@ -118,7 +118,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps enddo - B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 2d0*lambda*chi + B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi end do end do @@ -142,7 +142,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps enddo - B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 2d0*lambda*chi + B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi end do end do diff --git a/src/Makefile b/src/Makefile index 31d8569..7db1faa 100644 --- a/src/Makefile +++ b/src/Makefile @@ -52,7 +52,7 @@ debug: .DEFAULT_GOAL := default .PHONY: default debug -.PRECIOUS: $(LDIR)/%.a %/Makefile +.PRECIOUS: $(wildcard $(LDIR)/*.a) $(wildcard $(SDIR)/*/Makefile) clean: rm -f -- $(patsubst %/,$(LDIR)/%.a,$(wildcard */)) ; \ diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 29e771e..27206bf 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -356,7 +356,8 @@ program QuAcK ! Memory allocation - allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas)) + allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas), & + ERI_MO_bbbb(nBas,nBas,nBas,nBas),ERI_MO_abab(nBas,nBas,nBas,nBas)) ! 4-index transform for (aa|aa) block @@ -382,21 +383,14 @@ program QuAcK ket2 = 2 call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb) - if(spin_flip) then + ! 4-index transform for (ab|ab) block - allocate(ERI_MO_abab(nBas,nBas,nBas,nBas)) + bra1 = 1 + bra2 = 2 + ket1 = 1 + ket2 = 2 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_abab) - ! 4-index transform for (ab|ab) block - - bra1 = 1 - bra2 = 2 - ket1 = 1 - ket2 = 2 - call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_abab) - - end if - - else ! Memory allocation @@ -618,7 +612,7 @@ program QuAcK if(unrestricted) then call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_MO_aaaa,ERI_MO_aabb, & - ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) + ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF,cHF,S) else @@ -675,7 +669,7 @@ program QuAcK if(unrestricted) then call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF,cHF,S) else @@ -700,7 +694,7 @@ program QuAcK if(unrestricted) then call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF,cHF,S) else diff --git a/src/RPA/URPAx.f90 b/src/RPA/URPAx.f90 index 2f21783..33da34a 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/URPAx.f90 @@ -1,5 +1,5 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e,c,S) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -24,6 +24,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) @@ -87,8 +89,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) - call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - Omega_sc,XpY_sc,XmY_sc) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + c,S,Omega_sc,XpY_sc,XmY_sc) deallocate(Omega_sc,XpY_sc,XmY_sc) @@ -111,8 +113,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPAx ',6,nS_sf,Omega_sf) - call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & - Omega_sf,XpY_sf,XmY_sf) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & + c,S,Omega_sf,XpY_sf,XmY_sf) deallocate(Omega_sf,XpY_sf,XmY_sf) diff --git a/src/RPA/UdRPA.f90 b/src/RPA/UdRPA.f90 index 99d94de..1becf23 100644 --- a/src/RPA/UdRPA.f90 +++ b/src/RPA/UdRPA.f90 @@ -1,5 +1,5 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e,c,S) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -24,6 +24,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) @@ -86,8 +88,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPA ',5,nS_sc,Omega_sc) - call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - Omega_sc,XpY_sc,XmY_sc) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + c,S,Omega_sc,XpY_sc,XmY_sc) endif @@ -109,8 +111,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPA ',6,nS_sf,Omega_sf) - call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & - Omega_sf,XpY_sf,XmY_sf) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & + c,S,Omega_sf,XpY_sf,XmY_sf) endif diff --git a/src/utils/print_excitation.f90 b/src/utils/print_excitation.f90 new file mode 100644 index 0000000..c2bbdb4 --- /dev/null +++ b/src/utils/print_excitation.f90 @@ -0,0 +1,45 @@ +subroutine print_excitation(method,ispin,nS,Omega) + +! Print excitation energies for a given spin manifold + + implicit none + include 'parameters.h' + +! Input variables + + character*12,intent(in) :: method + integer,intent(in) :: ispin,nS + double precision,intent(in) :: Omega(nS) + +! Local variables + + character*14 :: spin_manifold + integer,parameter :: maxS = 32 + integer :: ia + + if(ispin == 1) spin_manifold = 'singlet' + if(ispin == 2) spin_manifold = 'triplet' + if(ispin == 3) spin_manifold = 'alpha-beta' + if(ispin == 4) spin_manifold = 'alpha-alpha' + if(ispin == 5) spin_manifold = 'spin-conserved' + if(ispin == 6) spin_manifold = 'spin-flip' + + write(*,*) + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + + do ia=1,min(nS,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------' + write(*,*) + +end subroutine print_excitation + +