Merge pull request #1155 from edoapra/bse-pull

BSE code
This commit is contained in:
NWChem: Open Source High-Performance Computational Chemistry 2025-10-17 22:17:21 -07:00 committed by GitHub
commit b5a513d7b0
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
64 changed files with 8568 additions and 1074 deletions

View file

@ -209,12 +209,6 @@ jobs:
armci_network: MPI-PR
nwchem_modules: "tinyqmpw python"
fc: ifx
- os: ubuntu-latest
experimental: true
mpi_impl: intel
armci_network: SOCKETS
nwchem_modules: "tce"
fc: ifx
- os: ubuntu-22.04
experimental: true
mpi_impl: intel

View file

@ -577,6 +577,9 @@ let "myexit+=$?"
let "myexit+=$?"
./runtests.mpi.unix procs $np gw_symmetry
let "myexit+=$?"
# BSE
./runtests.mpi.unix procs $np bse_ethene
let "myexit+=$?"
# x2c
./runtests.mpi.unix procs $np x2c-h2se
let "myexit+=$?"

View file

@ -433,6 +433,13 @@ sync
set overall_status = 1
continue
fi
# delete @GW lines for BSE tests
if [[ `grep -c 'NWChem BSE Module' ${STUB}.out` > 0 ]] ; then
rm -f sort.new
cat ${STUB}.ok.out.nwparse |grep -v @GW > sort.new && mv sort.new ${STUB}.ok.out.nwparse
rm -f sort.new
cat ${STUB}.out.nwparse |grep -v @GW > sort.new && mv sort.new ${STUB}.out.nwparse
fi
#check if output is from EOMCCSD, since EOMCCSD output is non-deterministic
if [[ `grep -c EOMCCSD ${STUB}.out.nwparse` > 0 ]] ; then
rm -f sort.new

View file

@ -0,0 +1,45 @@
title "ethene BSE"
memory total 2000 mb
start
set lindep:n_dep 0
geometry
H 0.000000 0.923274 1.238289
H 0.000000 -0.923274 1.238289
H 0.000000 0.923274 -1.238289
H 0.000000 -0.923274 -1.238289
C 0.000000 0.000000 0.668188
C 0.000000 0.000000 -0.668188
end
basis "ao basis" spherical bse
* library aug-cc-pvtz
end
basis "ri basis" spherical bse
* library aug-cc-pvtz-rifit
end
dft
direct
xc pbe0
tolerances acccoul 10
convergence energy 1d-8
grid fine
end
gw
evgw 5
method analytic
states alpha occ -1 vir -1
end
task gw
bse
nroots 12
method davidson
end
task bse

File diff suppressed because it is too large Load diff

View file

@ -35,11 +35,8 @@ end
gw
core
eta 0.01
method cdgw
states alpha occ 7
solver newton 15
convergence 0.001 ev
end
task dft gw
@ -52,27 +49,8 @@ end
gw
core
first 4
eta 0.01
method cdgw 128
solver newton 15
method cdgw
states alpha occ 4
convergence 0.001 ev
end
task dft gw
dft
noscf
end
unset gw:first
gw
core
eta 0.01
solver newton 15
states alpha occ 7
convergence 0.001 ev
end
task dft gw

View file

@ -1,4 +1,5 @@
argument 1 = gw_closedshell.nw
argument 1 = /people/meji656/calculations/bse/QA/gw_closedshell.nw
NWChem w/ OpenMP: maximum threads = 1
@ -21,11 +22,11 @@ H 2.5025 -1.2628 -0.0022
H 0.9511 -2.0241 -0.0008
end
basis "ao basis" spherical
basis "ao basis" bse spherical
* library cc-pvdz
end
basis "cd basis" spherical
basis "cd basis" bse spherical
* library cc-pvdz-rifit
end
@ -40,11 +41,8 @@ end
gw
core
eta 0.01
method cdgw
states alpha occ 7
solver newton 15
convergence 0.001 ev
end
task dft gw
@ -57,27 +55,8 @@ end
gw
core
first 4
eta 0.01
method cdgw 128
solver newton 15
method cdgw
states alpha occ 4
convergence 0.001 ev
end
task dft gw
dft
noscf
end
unset gw:first
gw
core
eta 0.01
solver newton 15
states alpha occ 7
convergence 0.001 ev
end
task dft gw
@ -88,7 +67,7 @@ task dft gw
Northwest Computational Chemistry Package (NWChem) 7.0.1
Northwest Computational Chemistry Package (NWChem) 7.3.0
--------------------------------------------------------
@ -96,7 +75,7 @@ task dft gw
Pacific Northwest National Laboratory
Richland, WA 99352
Copyright (c) 1994-2020
Copyright (c) 1994-2025
Pacific Northwest National Laboratory
Battelle Memorial Institute
@ -121,21 +100,21 @@ task dft gw
Job information
---------------
hostname = dmejiar-VirtualBox
program = nwchem
date = Fri Feb 4 04:13:35 2022
hostname = dc004
program = /qfs/people/meji656/software/nwchem.bse/bin/LINUX64/nwchem
date = Thu Oct 16 10:34:35 2025
compiled = Fri_Feb_04_04:04:58_2022
source = /home/dmejiar/Software/nwchem_gfortran_gw
nwchem branch = 7.0.0
nwchem revision = nwchem_on_git-3497-g0c7e480c37
ga revision = 5.8.1
use scalapack = T
input = gw_closedshell.nw
compiled = Thu_Oct_16_10:07:15_2025
source = /qfs/people/meji656/software/nwchem.bse
nwchem branch = 7.3.0
nwchem revision = v7.2.0-beta1-1624-g0dc011d0a0
ga revision = 5.9.2
use scalapack = F
input = /people/meji656/calculations/bse/QA/gw_closedshell.nw
prefix = gw_closedshell.
data base = ./gw_closedshell.db
status = startup
nproc = 2
nproc = 7
time left = -1s
@ -143,10 +122,10 @@ task dft gw
Memory information
------------------
heap = 65536000 doubles = 500.0 Mbytes
stack = 65535997 doubles = 500.0 Mbytes
global = 131072000 doubles = 1000.0 Mbytes (distinct from heap & stack)
total = 262143997 doubles = 2000.0 Mbytes
heap = 65535998 doubles = 500.0 Mbytes
stack = 65535995 doubles = 500.0 Mbytes
global = 131072000 doubles = 1000.0 Mbytes (distinct from heap & stack)
total = 262143993 doubles = 2000.0 Mbytes
verify = yes
hardfail = no
@ -311,9 +290,6 @@ task dft gw
library name resolved from: environment
library file name is: </home/dmejiar/Software/nwchem_gfortran_gw/src/basis/libraries.bse/>
Summary of "ao basis" -> "" (spherical)
@ -323,9 +299,6 @@ task dft gw
* cc-pvdz on all atoms
library name resolved from: environment
library file name is: </home/dmejiar/Software/nwchem_gfortran_gw/src/basis/libraries.bse/>
Summary of "cd basis" -> "" (spherical)
@ -337,8 +310,8 @@ task dft gw
pbe96 is a nonlocal functional; adding pw91lda local functional.
NWChem GW@DFT Module
--------------------
NWChem GW Module
----------------
GW Trifluoroacetamide
@ -796,48 +769,48 @@ task dft gw
HOMO = -0.354463
LUMO = 0.023799
Time after variat. SCF: 1.2
Time after variat. SCF: 0.8
3 Center 2 Electron Integral Information
----------------------------------------
Maximum number of 3-center 2e- integrals is: 4898880.
This is reduced with Schwarz screening to: 3084900.
Incore requires a per proc buffer size of: 1377181.
Incore requires a per proc buffer size of: 398581.
Minimum dble words available (all nodes) is: 131060988
This is reduced (for later use) to: 130785792
proc 0 Suggested buffer size is: 1377181
Max Suggested buffer size is: 1377181
This is reduced (for later use) to: 130844112
proc 0 Suggested buffer size is: 398581
Max Suggested buffer size is: 398581
no. integral batches is: 1000
Incore memory use for 3-center 2e- integrals is turned off.
Time prior to 1st pass: 1.2
Time prior to 1st pass: 0.8
Memory utilization after 1st SCF pass:
Heap Space remaining (MW): 65.53 65525792
Stack Space remaining (MW): 65.54 65535156
Heap Space remaining (MW): 65.53 65525790
Stack Space remaining (MW): 65.54 65535148
convergence iter energy DeltaE RMS-Dens Diis-err time
---------------- ----- ----------------- --------- --------- --------- ------
d= 0,ls=0.0,diis 1 -506.2138586363 -8.41D+02 1.85D-02 2.46D+00 5.4
d= 0,ls=0.0,diis 2 -505.9428729180 2.71D-01 1.15D-02 5.90D+00 9.7
d= 0,ls=0.0,diis 3 -506.4729357426 -5.30D-01 2.21D-03 1.04D-01 13.9
d= 0,ls=0.0,diis 4 -506.4807116964 -7.78D-03 8.40D-04 2.68D-02 18.1
d= 0,ls=0.0,diis 5 -506.4832631842 -2.55D-03 1.27D-04 5.68D-04 23.8
d= 0,ls=0.0,diis 1 -506.2138586373 -8.41D+02 1.85D-02 2.46D+00 1.3
d= 0,ls=0.0,diis 2 -505.9428729965 2.71D-01 1.15D-02 5.90D+00 1.9
d= 0,ls=0.0,diis 3 -506.4729357440 -5.30D-01 2.21D-03 1.04D-01 2.4
d= 0,ls=0.0,diis 4 -506.4807116867 -7.78D-03 8.40D-04 2.68D-02 3.0
d= 0,ls=0.0,diis 5 -506.4832631820 -2.55D-03 1.27D-04 5.68D-04 3.7
Resetting Diis
d= 0,ls=0.0,diis 6 -506.4833170187 -5.38D-05 4.21D-05 2.12D-05 29.7
d= 0,ls=0.0,diis 7 -506.4833201162 -3.10D-06 2.25D-05 3.93D-06 36.1
d= 0,ls=0.0,diis 8 -506.4833202140 -9.77D-08 9.51D-06 4.40D-06 42.1
d= 0,ls=0.0,diis 6 -506.4833170167 -5.38D-05 4.21D-05 2.12D-05 4.4
d= 0,ls=0.0,diis 7 -506.4833201141 -3.10D-06 2.25D-05 3.93D-06 5.2
d= 0,ls=0.0,diis 8 -506.4833202118 -9.77D-08 9.51D-06 4.40D-06 5.9
Total DFT energy = -506.483320213950
One electron energy = -1365.175909550851
Coulomb energy = 581.765579688889
Exchange-Corr. energy = -58.042890125742
Total DFT energy = -506.483320211761
One electron energy = -1365.175909547848
Coulomb energy = 581.765579681333
Exchange-Corr. energy = -58.042890118999
Nuclear repulsion energy = 334.969899773753
Numeric. integr. density = 56.000003719920
Numeric. integr. density = 56.000003716419
Total iterative time = 40.9s
Total iterative time = 5.1s
@ -870,32 +843,36 @@ task dft gw
2 0 0 2 -28.307350 -52.360374 -52.360374 76.413398
***********************
*** GW@DFT ***
***********************
***********************
*** GW ***
***********************
Memory Information
------------------
Available GA per process is 999.9 MBs
Available GA per process is 1000.0 MBs
Available stack per process is 500.0 MBs
Available heap per process is 499.9 MBs
Peak memory in Initialization
Peak GA per process is 5.6 MBs
Peak GA per process is 1.6 MBs
Peak stack per process is 0.1 MBs
Peak heap per process is 6.5 MBs
Peak heap per process is 2.9 MBs
Memory not in MA per process 0.0 MBs
Peak memory in CDGW Calculation
Peak GA per process is 7.0 MBs
Peak GA per process is 2.0 MBs
Peak stack per process is 2.7 MBs
Peak heap per process is 2.1 MBs
Peak heap per process is 1.3 MBs
Memory not in MA per process 0.0 MBs
Computing 2-center integrals
Computing Inverse Cholesky factor
Computing 3-center integrals
Transform ERIs
Distribute ERIs
Read inverse
Orthogonalize ERIs
Computing Vxc
Computing Sigma_x
Computing W(iw) on the imaginary grid
@ -904,28 +881,28 @@ task dft gw
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 1 -695.048 0.000
@GW 2 -694.653 0.000
@GW 3 -694.648 0.000
@GW 4 -538.987 0.000
@GW 5 -408.010 0.000
@GW 6 -300.254 0.000
@GW 7 -296.783 0.000
@GW 1 -695.046 0.002
@GW 2 -694.650 0.002
@GW 3 -694.645 0.002
@GW 4 -538.980 0.006
@GW 5 -408.006 0.001
@GW 6 -300.252 0.002
@GW 7 -296.781 0.002
GW Timing Statistics
------------------------
Init. : 0.1s
ERIs : 0.7s
ERIs : 0.2s
Distr.: 0.0s
Vxc : 4.9s
Vxc : 0.5s
Vmo : 0.0s
iW : 2.8s
iW : 0.9s
I_n : 0.0s
R_n : 1.1s
GW total: 9.6s
R_n : 0.2s
GW total: 2.0s
Task times cpu: 51.4s wall: 51.6s
Task times cpu: 7.6s wall: 7.6s
NWChem Input Module
@ -933,8 +910,8 @@ task dft gw
NWChem GW@DFT Module
--------------------
NWChem GW Module
----------------
GW Trifluoroacetamide
@ -1034,34 +1011,34 @@ task dft gw
Loading old vectors from job with title :
GW Trifluoroacetamide
W vals
Time after variat. SCF: 51.9
Time after variat. SCF: 8.0
3 Center 2 Electron Integral Information
----------------------------------------
Maximum number of 3-center 2e- integrals is: 4898880.
This is reduced with Schwarz screening to: 3084900.
Incore requires a per proc buffer size of: 1377181.
Incore requires a per proc buffer size of: 398581.
Minimum dble words available (all nodes) is: 131060988
This is reduced (for later use) to: 130586152
proc 0 Suggested buffer size is: 1377181
Max Suggested buffer size is: 1377181
This is reduced (for later use) to: 130644472
proc 0 Suggested buffer size is: 398581
Max Suggested buffer size is: 398581
no. integral batches is: 1000
Incore memory use for 3-center 2e- integrals is turned off.
Time prior to 1st pass: 51.9
Time prior to 1st pass: 8.0
Total DFT energy = -506.483320618786
One electron energy = -1365.176616645610
Coulomb energy = 581.766179275143
Exchange-Corr. energy = -58.042783022072
Total DFT energy = -506.483320616637
One electron energy = -1365.176616647789
Coulomb energy = 581.766179277314
Exchange-Corr. energy = -58.042783019915
Nuclear repulsion energy = 334.969899773753
Numeric. integr. density = 56.000003720190
Numeric. integr. density = 56.000003716815
Total iterative time = 5.9s
Total iterative time = 0.7s
@ -1094,32 +1071,36 @@ GW Trifluoroacetamide
2 0 0 2 -28.307350 -52.360374 -52.360374 76.413398
***********************
*** GW@DFT ***
***********************
***********************
*** GW ***
***********************
Memory Information
------------------
Available GA per process is 999.9 MBs
Available GA per process is 1000.0 MBs
Available stack per process is 500.0 MBs
Available heap per process is 499.9 MBs
Peak memory in Initialization
Peak GA per process is 5.6 MBs
Peak GA per process is 1.6 MBs
Peak stack per process is 0.1 MBs
Peak heap per process is 6.5 MBs
Peak heap per process is 2.9 MBs
Memory not in MA per process 0.0 MBs
Peak memory in CDGW Calculation
Peak GA per process is 7.0 MBs
Peak GA per process is 2.0 MBs
Peak stack per process is 2.7 MBs
Peak heap per process is 1.5 MBs
Peak heap per process is 1.2 MBs
Memory not in MA per process 0.0 MBs
Computing 2-center integrals
Computing Inverse Cholesky factor
Computing 3-center integrals
Transform ERIs
Distribute ERIs
Read inverse
Orthogonalize ERIs
Computing Vxc
Computing Sigma_x
Computing W(iw) on the imaginary grid
@ -1128,249 +1109,25 @@ GW Trifluoroacetamide
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 4 -538.987 0.000
@GW 5 -408.010 0.000
@GW 6 -300.254 0.000
@GW 7 -296.783 0.000
@GW 4 -538.980 0.006
@GW 5 -408.006 0.001
@GW 6 -300.252 0.002
@GW 7 -296.781 0.002
GW Timing Statistics
------------------------
Init. : 0.1s
ERIs : 0.8s
ERIs : 0.2s
Distr.: 0.0s
Vxc : 4.9s
Vxc : 0.5s
Vmo : 0.0s
iW : 1.6s
iW : 0.9s
I_n : 0.0s
R_n : 0.6s
GW total: 8.0s
R_n : 0.1s
GW total: 1.8s
Task times cpu: 14.1s wall: 14.1s
NWChem Input Module
-------------------
NWChem GW@DFT Module
--------------------
GW Trifluoroacetamide
NWChem DFT Module
-----------------
GW Trifluoroacetamide
Summary of "ao basis" -> "ao basis" (spherical)
------------------------------------------------------------------------------
Tag Description Shells Functions and Types
---------------- ------------------------------ ------ ---------------------
F cc-pvdz 6 14 3s2p1d
O cc-pvdz 6 14 3s2p1d
N cc-pvdz 6 14 3s2p1d
C cc-pvdz 6 14 3s2p1d
H cc-pvdz 3 5 2s1p
Caching 1-el integrals
General Information
-------------------
SCF calculation type: DFT
Wavefunction type: closed shell.
No. of atoms : 9
No. of electrons : 56
Alpha electrons : 28
Beta electrons : 28
Charge : 0
Spin multiplicity: 1
Use of symmetry is: off; symmetry adaption is: off
Maximum number of iterations: 50
This is a Direct SCF calculation.
AO basis - number of functions: 108
number of shells: 48
A Charge density fitting basis will be used.
CD basis - number of functions: 420
number of shells: 138
Convergence on energy requested: 1.00D-06
Convergence on density requested: 1.00D-05
Convergence on gradient requested: 5.00D-04
XC Information
--------------
Hartree-Fock (Exact) Exchange 0.450
PerdewBurkeErnzerhof Exchange Functional 0.550
Perdew 1991 LDA Correlation Functional 1.000 local
PerdewBurkeErnz. Correlation Functional 1.000 non-local
Grid Information
----------------
Grid used for XC integration: medium
Radial quadrature: Mura-Knowles
Angular quadrature: Lebedev.
Tag B.-S. Rad. Rad. Pts. Rad. Cut. Ang. Pts.
--- ---------- --------- --------- ---------
F 0.50 49 11.0 434
O 0.60 49 9.0 434
N 0.65 49 8.0 434
C 0.70 49 10.0 434
H 0.35 45 10.0 434
Grid pruning is: on
Number of quadrature shells: 433
Spatial weights used: Erf1
Convergence Information
-----------------------
Convergence aids based upon iterative change in
total energy or number of iterations.
Levelshifting, if invoked, occurs when the
HOMO/LUMO gap drops below (HL_TOL): 1.00D-02
DIIS, if invoked, will attempt to extrapolate
using up to (NFOCK): 10 stored Fock matrices.
Damping( 0%) Levelshifting(0.5) DIIS
--------------- ------------------- ---------------
dE on: start ASAP start
dE off: 2 iters 50 iters 50 iters
Screening Tolerance Information
-------------------------------
Density screening/tol_rho: 1.00D-10
AO Gaussian exp screening on grid/accAOfunc: 14
CD Gaussian exp screening on grid/accCDfunc: 20
XC Gaussian exp screening on grid/accXCfunc: 20
Schwarz screening/accCoul: 1.00D-12
Loading old vectors from job with title :
GW Trifluoroacetamide
Time after variat. SCF: 66.0
3 Center 2 Electron Integral Information
----------------------------------------
Maximum number of 3-center 2e- integrals is: 4898880.
This is reduced with Schwarz screening to: 3084900.
Incore requires a per proc buffer size of: 1377181.
Minimum dble words available (all nodes) is: 131060988
This is reduced (for later use) to: 130586152
proc 0 Suggested buffer size is: 1377181
Max Suggested buffer size is: 1377181
no. integral batches is: 1000
Incore memory use for 3-center 2e- integrals is turned off.
Time prior to 1st pass: 66.0
Total DFT energy = -506.483320619392
One electron energy = -1365.176616645610
Coulomb energy = 581.766179275142
Exchange-Corr. energy = -58.042783022677
Nuclear repulsion energy = 334.969899773753
Numeric. integr. density = 56.000003721840
Total iterative time = 5.9s
center of mass
--------------
x = 0.05046829 y = 0.11377991 z = 0.00004386
moments of inertia (a.u.)
------------------
724.808949354496 0.000000000000 0.000000000000
0.000000000000 484.336918953950 -0.000000000000
0.000000000000 -0.000000000000 886.538753836127
Multipole analysis of the density
---------------------------------
L x y z total alpha beta nuclear
- - - - ----- ----- ---- -------
0 0 0 0 -0.000000 -28.000000 -28.000000 56.000000
1 1 0 0 -1.489094 -0.744547 -0.744547 -0.000000
1 0 1 0 -0.546971 -0.273486 -0.273486 0.000000
1 0 0 1 0.000605 0.000302 0.000302 0.000000
2 2 0 0 -27.181301 -96.832025 -96.832025 166.482750
2 1 1 0 4.886093 -1.980672 -1.980672 8.847437
2 1 0 1 0.002838 -0.001223 -0.001223 0.005284
2 0 2 0 -26.258758 -158.160655 -158.160655 290.062553
2 0 1 1 -0.001360 -0.006340 -0.006340 0.011320
2 0 0 2 -28.307350 -52.360374 -52.360374 76.413398
***********************
*** GW@DFT ***
***********************
Memory Information
------------------
Available GA per process is 999.9 MBs
Available stack per process is 500.0 MBs
Available heap per process is 499.9 MBs
Peak memory in Initialization
Peak GA per process is 5.6 MBs
Peak stack per process is 0.1 MBs
Peak heap per process is 6.5 MBs
Memory not in MA per process 0.0 MBs
Peak memory in AFF Calculation
Peak GA per process is 27.6 MBs
Peak stack per process is 38.3 MBs
Peak heap per process is 0.0 MBs
Memory not in MA per process 38.3 MBs
Computing 2-center integrals
Computing Inverse Cholesky factor
Computing 3-center integrals
Computing V_xc ...
Computing Sigma_x
Computing RPA polarizability
Solving quasiparticle equations
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 1 -695.048 0.000
@GW 2 -694.652 0.000
@GW 3 -694.648 0.000
@GW 4 -538.986 0.000
@GW 5 -408.009 0.000
@GW 6 -300.254 0.000
@GW 7 -296.783 0.000
GW Timing Statistics
------------------------
Init. : 0.1s
ERIs : 0.7s
Distr.: 0.0s
Vxc : 4.7s
Sigma_x : 0.0s
RPA : 2.6s
Wmn : 0.1s
Sigma_c : 0.0s
GW total: 8.4s
Task times cpu: 14.3s wall: 14.4s
Task times cpu: 2.6s wall: 2.6s
NWChem Input Module
@ -1390,11 +1147,11 @@ MA usage statistics:
heap stack
---- -----
current number of blocks 0 0
maximum number of blocks 28 56
maximum number of blocks 28 62
current total bytes 0 0
maximum total bytes 25980280 80584200
maximum total K-bytes 25981 80585
maximum total M-bytes 26 81
maximum total bytes 22153480 30740728
maximum total K-bytes 22154 30741
maximum total M-bytes 23 31
CITATION
@ -1430,24 +1187,25 @@ MA usage statistics:
AUTHORS
-------
E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, M. Valiev, W. A. de Jong,
T. P. Straatsma, H. J. J. van Dam, D. Wang, T. L. Windus, N. P. Bauman,
A. Panyala, J. Hammond, J. Autschbach, K. Bhaskaran-Nair, J. Brabec,
K. Lopata, S. A. Fischer, S. Krishnamoorthy, M. Jacquelin, W. Ma, M. Klemm,
O. Villa, Y. Chen, V. Anisimov, F. Aquino, S. Hirata, M. T. Hackler,
Eric Hermes, L. Jensen, J. E. Moore, J. C. Becca, V. Konjkov,
D. Mejia-Rodriguez, T. Risthaus, M. Malagoli, A. Marenich,
A. Otero-de-la-Roza, J. Mullin, P. Nichols, R. Peverati, J. Pittner, Y. Zhao,
P.-D. Fan, A. Fonari, M. J. Williamson, R. J. Harrison, J. R. Rehr,
M. Dupuis, D. Silverstein, D. M. A. Smith, J. Nieplocha, V. Tipparaju,
M. Krishnan, B. E. Van Kuiken, A. Vazquez-Mayagoitia, M. Swart, Q. Wu,
T. Van Voorhis, A. A. Auer, M. Nooijen, L. D. Crosby, E. Brown, G. Cisneros,
G. I. Fann, H. Fruchtl, J. Garza, K. Hirao, R. A. Kendall, J. A. Nichols,
K. Tsemekhman, K. Wolinski, J. Anchell, D. E. Bernholdt, P. Borowski,
T. Clark, D. Clerc, H. Dachsel, M. J. O. Deegan, K. Dyall, D. Elwood,
E. Glendening, M. Gutowski, A. C. Hess, J. Jaffe, B. G. Johnson, J. Ju,
R. Kobayashi, R. Kutteh, Z. Lin, R. Littlefield, X. Long, B. Meng,
T. Nakajima, S. Niu, L. Pollack, M. Rosing, K. Glaesemann, G. Sandrone,
M. Stave, H. Taylor, G. Thomas, J. H. van Lenthe, A. T. Wong, Z. Zhang.
E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, D. Mejia-Rodriguez,
N. P. Bauman, A. Panyala, R. J. Harrison, M. Valiev, W. A. de Jong,
T. P. Straatsma, H. J. J. van Dam, D. Wang, T. L. Windus, J. Hammond,
J. Autschbach, R. de P. Soares, A. Kunitsa, K. Bhaskaran-Nair, J. Brabec,
K. Lopata, S. A. Fischer, S. Krishnamoorthy, M. Jacquelin, R. Goswami,
A. Woods, W. Ma, M. Klemm, O. Villa, Y. Chen, V. Anisimov, F. Aquino,
S. Hirata, M. T. Hackler, E. Hermes, L. Jensen, J. E. Moore, J. C. Becca,
V. Konjkov, T. Risthaus, M. Malagoli, A. Marenich, A. Otero-de-la-Roza,
J. Mullin, P. Nichols, R. Peverati, J. Pittner, Y. Zhao, P.-D. Fan,
A. Fonari, M. J. Williamson, J. R. Rehr, M. Dupuis, D. Silverstein,
D. M. A. Smith, J. Nieplocha, V. Tipparaju, M. Krishnan, B. E. Van Kuiken,
A. Vazquez-Mayagoitia, M. Swart, Q. Wu, T. Van Voorhis, A. A. Auer,
M. Nooijen, L. D. Crosby, E. Brown, G. Cisneros, G. I. Fann, H. Fruchtl,
J. Garza, K. Hirao, R. A. Kendall, J. A. Nichols, K. Tsemekhman, K. Wolinski,
J. Anchell, D. E. Bernholdt, P. Borowski, T. Clark, D. Clerc, H. Dachsel,
M. J. O. Deegan, K. Dyall, D. Elwood, E. Glendening, M. Gutowski, A. C. Hess,
J. Jaffe, B. G. Johnson, J. Ju, R. Kobayashi, R. Kutteh, Z. Lin,
R. Littlefield, X. Long, B. Meng, T. Nakajima, S. Niu, L. Pollack, M. Rosing,
K. Glaesemann, G. Sandrone, M. Stave, H. Taylor, G. Thomas, J. H. van Lenthe,
A. T. Wong, Z. Zhang.
Total times cpu: 79.8s wall: 80.2s
Total times cpu: 10.5s wall: 10.5s

View file

@ -24,35 +24,18 @@ end
dft
odft
mult 2
xc pbe96
xc pbe0
direct
grid fine nodisk
tolerances accCoul 14
tolerances accCoul 12
convergence energy 1d-8
noprint "final vectors analysis"
end
gw
eta 0.001
method cdgw
solver newton 20
states alpha occ 5 vir 5
states beta occ 5 vir 5
convergence 0.001 ev
end
task dft gw
dft
noscf
end
gw
eta 0.001
solver newton 20
states alpha occ 5 vir 5
states beta occ 5 vir 5
convergence 0.001 ev
end
task dft gw

View file

@ -1,4 +1,5 @@
argument 1 = gw_openshell.nw
argument 1 = /data/edo/nwchem/nwchem_bse-pull//QA/tests/gw_openshell/gw_openshell.nw
NWChem w/ OpenMP: maximum threads = 1
@ -17,11 +18,11 @@ H 2.5025 -1.2628 -0.0022
H 0.9511 -2.0241 -0.0008
end
basis "ao basis" spherical
basis "ao basis" bse spherical
* library cc-pvdz
end
basis "cd basis" spherical
basis "cd basis" bse spherical
* library cc-pvdz-rifit
end
@ -29,35 +30,18 @@ end
dft
odft
mult 2
xc pbe96
xc pbe0
direct
grid fine nodisk
tolerances accCoul 14
tolerances accCoul 12
convergence energy 1d-8
noprint "final vectors analysis"
end
gw
eta 0.001
method cdgw
solver newton 20
states alpha occ 5 vir 5
states beta occ 5 vir 5
convergence 0.001 ev
end
task dft gw
dft
noscf
end
gw
eta 0.001
solver newton 20
states alpha occ 5 vir 5
states beta occ 5 vir 5
convergence 0.001 ev
end
task dft gw
@ -69,7 +53,7 @@ task dft gw
Northwest Computational Chemistry Package (NWChem) 7.0.1
Northwest Computational Chemistry Package (NWChem) 7.3.0
--------------------------------------------------------
@ -77,7 +61,7 @@ task dft gw
Pacific Northwest National Laboratory
Richland, WA 99352
Copyright (c) 1994-2020
Copyright (c) 1994-2025
Pacific Northwest National Laboratory
Battelle Memorial Institute
@ -102,21 +86,21 @@ task dft gw
Job information
---------------
hostname = dmejiar-VirtualBox
program = nwchem
date = Fri Feb 4 10:51:34 2022
hostname = mica
program = /data/edo/nwchem/nwchem_bse-pull//bin/LINUX64/nwchem
date = Thu Oct 16 12:25:05 2025
compiled = Fri_Feb_04_04:17:52_2022
source = /home/dmejiar/Software/nwchem_gfortran_gw
nwchem branch = 7.0.0
nwchem revision = nwchem_on_git-3497-g0c7e480c37
ga revision = 5.8.1
compiled = Thu_Oct_16_12:11:32_2025
source = /data/edo/nwchem/nwchem_bse-pull/
nwchem branch = 7.3.0
nwchem revision = v7.2.0-beta1-1593-g3a8d752956
ga revision = 5.9.2
use scalapack = T
input = gw_openshell.nw
input = /data/edo/nwchem/nwchem_bse-pull//QA/tests/gw_openshell/gw_openshell.nw
prefix = gw_openshell.
data base = ./gw_openshell.db
status = startup
nproc = 2
nproc = 1
time left = -1s
@ -124,10 +108,10 @@ task dft gw
Memory information
------------------
heap = 65535996 doubles = 500.0 Mbytes
stack = 65536001 doubles = 500.0 Mbytes
global = 131072000 doubles = 1000.0 Mbytes (distinct from heap & stack)
total = 262143997 doubles = 2000.0 Mbytes
heap = 65535996 doubles = 500.0 Mbytes
stack = 65536001 doubles = 500.0 Mbytes
global = 131072000 doubles = 1000.0 Mbytes (distinct from heap & stack)
total = 262143997 doubles = 2000.0 Mbytes
verify = yes
hardfail = no
@ -251,9 +235,6 @@ task dft gw
library name resolved from: environment
library file name is: </home/dmejiar/Software/nwchem_gfortran_gw/src/basis/libraries.bse/>
Summary of "ao basis" -> "" (spherical)
@ -263,9 +244,6 @@ task dft gw
* cc-pvdz on all atoms
library name resolved from: environment
library file name is: </home/dmejiar/Software/nwchem_gfortran_gw/src/basis/libraries.bse/>
Summary of "cd basis" -> "" (spherical)
@ -276,8 +254,8 @@ task dft gw
NWChem GW@DFT Module
--------------------
NWChem GW Module
----------------
GW Trifluoroacetamide
@ -603,8 +581,9 @@ task dft gw
XC Information
--------------
PBE96 Method XC Functional
PerdewBurkeErnzerhof Exchange Functional 1.000
PBE0 Method XC Functional
Hartree-Fock (Exact) Exchange 0.250
PerdewBurkeErnzerhof Exchange Functional 0.750
Perdew 1991 LDA Correlation Functional 1.000 local
PerdewBurkeErnz. Correlation Functional 1.000 non-local
@ -644,7 +623,7 @@ task dft gw
AO Gaussian exp screening on grid/accAOfunc: 18
CD Gaussian exp screening on grid/accCDfunc: 20
XC Gaussian exp screening on grid/accXCfunc: 20
Schwarz screening/accCoul: 1.00D-14
Schwarz screening/accCoul: 1.00D-12
Superposition of Atomic Density Guess
@ -661,21 +640,21 @@ task dft gw
HOMO = -0.154097
LUMO = 0.070838
Time after variat. SCF: 0.6
Time after variat. SCF: 0.7
3 Center 2 Electron Integral Information
----------------------------------------
Maximum number of 3-center 2e- integrals is: 529984.
This is reduced with Schwarz screening to: 367500.
Incore requires a per proc buffer size of: 175029.
Incore requires a per proc buffer size of: 327517.
Minimum dble words available (all nodes) is: 131065476
This is reduced (for later use) to: 130813264
proc 0 Suggested buffer size is: 175029
Max Suggested buffer size is: 175029
This is reduced (for later use) to: 130794336
proc 0 Suggested buffer size is: 327517
Max Suggested buffer size is: 327517
no. integral batches is: 1000
Incore memory use for 3-center 2e- integrals is turned off.
Time prior to 1st pass: 0.6
Time prior to 1st pass: 0.7
Memory utilization after 1st SCF pass:
Heap Space remaining (MW): 65.53 65530108
@ -683,42 +662,42 @@ task dft gw
convergence iter energy DeltaE RMS-Dens Diis-err time
---------------- ----- ----------------- --------- --------- --------- ------
d= 0,ls=0.0,diis 1 -168.9561121185 -2.32D+02 2.05D-02 2.18D-01 1.5
1.84D-02 2.17D-01
d= 0,ls=0.0,diis 2 -167.9504224824 1.01D+00 1.76D-02 2.63D+00 2.5
1.47D-02 2.33D+00
d= 0,ls=0.0,diis 3 -169.0537934402 -1.10D+00 3.84D-03 2.91D-02 3.3
2.75D-03 2.86D-02
d= 0,ls=0.0,diis 4 -169.0564442937 -2.65D-03 2.06D-03 2.17D-02 4.1
1.46D-03 2.00D-02
d= 0,ls=0.0,diis 5 -169.0660484759 -9.60D-03 2.88D-04 5.93D-04 4.8
2.45D-04 5.80D-04
d= 0,ls=0.0,diis 1 -168.9657539780 -2.32D+02 1.39D-02 1.93D-01 1.6
1.31D-02 1.92D-01
d= 0,ls=0.0,diis 2 -168.6178261291 3.48D-01 1.00D-02 1.13D+00 2.5
9.06D-03 1.06D+00
d= 0,ls=0.0,diis 3 -169.0613313213 -4.44D-01 1.36D-03 7.44D-03 3.5
1.32D-03 8.46D-03
d= 0,ls=0.0,diis 4 -169.0630054782 -1.67D-03 9.29D-04 4.19D-03 4.5
5.76D-04 3.46D-03
d= 0,ls=0.0,diis 5 -169.0646959497 -1.69D-03 1.55D-04 1.38D-04 5.5
1.31D-04 1.63D-04
Resetting Diis
d= 0,ls=0.0,diis 6 -169.0663002306 -2.52D-04 1.01D-04 2.23D-05 5.6
6.70D-05 2.22D-05
d= 0,ls=0.0,diis 7 -169.0663113153 -1.11D-05 8.27D-05 2.77D-06 6.5
5.41D-05 1.73D-06
d= 0,ls=0.0,diis 8 -169.0663035508 7.76D-06 4.90D-05 2.02D-05 7.2
4.55D-05 2.15D-05
d= 0,ls=0.0,diis 9 -169.0663123641 -8.81D-06 1.02D-05 3.98D-07 8.1
7.29D-06 3.32D-07
d= 0,ls=0.0,diis 10 -169.0663125308 -1.67D-07 2.91D-06 3.07D-08 8.9
2.57D-06 3.01D-08
d= 0,ls=0.0,diis 11 -169.0663125423 -1.15D-08 1.15D-06 1.20D-08 9.6
1.36D-06 1.15D-08
d= 0,ls=0.0,diis 12 -169.0663125471 -4.87D-09 6.52D-07 1.79D-09 10.4
2.64D-07 1.06D-09
d= 0,ls=0.0,diis 6 -169.0647672024 -7.13D-05 4.96D-05 7.61D-06 6.5
5.91D-05 7.90D-06
d= 0,ls=0.0,diis 7 -169.0647730165 -5.81D-06 5.03D-05 2.16D-06 7.5
1.81D-05 8.25D-07
d= 0,ls=0.0,diis 8 -169.0647729683 4.81D-08 1.78D-05 2.36D-06 8.4
2.14D-05 3.64D-06
d= 0,ls=0.0,diis 9 -169.0647741451 -1.18D-06 9.22D-06 5.57D-07 9.4
4.55D-06 2.30D-07
d= 0,ls=0.0,diis 10 -169.0647743526 -2.08D-07 9.00D-06 4.08D-08 10.4
8.19D-06 2.42D-08
d= 0,ls=0.0,diis 11 -169.0647744288 -7.62D-08 2.13D-06 7.37D-09 11.4
1.35D-06 1.44D-09
d= 0,ls=0.0,diis 12 -169.0647744339 -5.12D-09 7.95D-07 7.61D-10 12.4
6.20D-07 8.83D-10
Total DFT energy = -169.066312547124
One electron energy = -355.393819271438
Coulomb energy = 144.424914623888
Exchange-Corr. energy = -21.405786398232
Total DFT energy = -169.064774433938
One electron energy = -355.358774140560
Coulomb energy = 144.394070505297
Exchange-Corr. energy = -21.408449297333
Nuclear repulsion energy = 63.308378498658
Numeric. integr. density = 22.999999954116
Numeric. integr. density = 22.999999950854
Total iterative time = 9.8s
Total iterative time = 11.7s
@ -728,27 +707,27 @@ task dft gw
alpha 1 2 3 4 5 6 7 8 9 10
beta 1 2 3 4 5 6 7 8 9 10
overlap 1.000 1.000 1.000 1.000 1.000 0.996 0.996 0.990 1.000 0.985
overlap 1.000 1.000 1.000 1.000 0.999 0.993 0.991 0.974 1.000 0.968
alpha 11 12 13 14 15 16 17 18 19 20
beta 11 12 13 14 15 16 17 18 19 20
overlap 1.000 0.993 0.999 0.994 0.995 0.999 0.999 0.952 1.000 0.945
overlap 1.000 0.985 0.999 0.987 0.993 0.998 0.997 0.946 1.000 0.920
alpha 21 22 23 24 25 26 27 28 29 30
beta 21 22 23 24 25 26 27 28 29 30
overlap 0.958 0.999 0.998 1.000 0.998 0.999 0.998 0.996 0.997 0.980
beta 21 22 23 24 25 26 27 28 30 29
overlap 0.928 0.997 0.998 1.000 0.997 0.999 0.996 0.993 0.996 0.977
alpha 31 32 33 34 35 36 37 38 39 40
beta 31 32 33 34 35 36 37 38 39 40
overlap 0.985 0.996 0.991 0.992 0.999 1.000 0.997 0.994 0.997 0.999
overlap 0.983 0.995 0.992 0.993 0.999 1.000 0.986 0.979 0.993 1.000
alpha 41 42 43 44 45 46 47 48 49 50
beta 41 42 43 44 45 46 47 48 49 50
overlap 0.999 0.996 0.994 0.999 0.998 0.998 1.000 0.999 0.999 1.000
overlap 0.999 0.999 0.992 0.994 0.998 0.997 1.000 0.998 0.999 1.000
alpha 51 52
@ -758,7 +737,7 @@ task dft gw
--------------------------
Expectation value of S2:
--------------------------
<S2> = 0.7519 (Exact = 0.7500)
<S2> = 0.7537 (Exact = 0.7500)
center of mass
@ -778,21 +757,21 @@ task dft gw
- - - - ----- ----- ---- -------
0 0 0 0 -0.000000 -12.000000 -11.000000 23.000000
1 1 0 0 0.238427 0.541429 -0.303002 0.000000
1 0 1 0 -1.321737 -0.969708 -0.352030 0.000000
1 0 0 1 0.000663 -0.000054 0.000717 0.000000
1 1 0 0 0.206401 0.547712 -0.341311 0.000000
1 0 1 0 -1.390130 -1.013204 -0.376927 0.000000
1 0 0 1 0.000816 0.000001 0.000815 0.000000
2 2 0 0 -11.884337 -13.068472 -9.699498 10.883633
2 1 1 0 -0.250985 0.506239 0.213062 -0.970286
2 1 0 1 0.006019 -0.000243 -0.000755 0.007017
2 0 2 0 -11.275818 -52.448505 -49.703057 90.875744
2 0 1 1 0.001863 0.000480 -0.000509 0.001891
2 0 0 2 -13.100006 -6.900236 -6.199822 0.000052
2 2 0 0 -11.728128 -12.959803 -9.651957 10.883633
2 1 1 0 -0.304140 0.447969 0.218176 -0.970286
2 1 0 1 0.006364 0.000093 -0.000746 0.007017
2 0 2 0 -11.396501 -52.445939 -49.826305 90.875744
2 0 1 1 0.001878 0.000424 -0.000437 0.001891
2 0 0 2 -13.115703 -6.888717 -6.227038 0.000052
***********************
*** GW@DFT ***
***********************
***********************
*** GW ***
***********************
Memory Information
@ -802,15 +781,15 @@ task dft gw
Available heap per process is 500.0 MBs
Peak memory in Initialization
Peak GA per process is 1.4 MBs
Peak GA per process is 2.7 MBs
Peak stack per process is 0.0 MBs
Peak heap per process is 1.6 MBs
Peak heap per process is 2.9 MBs
Memory not in MA per process 0.0 MBs
Peak memory in CDGW Calculation
Peak GA per process is 1.7 MBs
Peak GA per process is 3.4 MBs
Peak stack per process is 0.6 MBs
Peak heap per process is 1.5 MBs
Peak heap per process is 1.9 MBs
Memory not in MA per process 0.0 MBs
Computing 2-center integrals
@ -824,323 +803,48 @@ task dft gw
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 8 -13.232 0.003
@GW 9 -12.555 0.001
@GW 10 -12.357 0.000
@GW 11 -9.991 0.000
@GW 12 -7.779 0.001
@GW 13 3.732 0.000
@GW 14 3.851 0.001
@GW 15 5.314 0.000
@GW 16 9.836 0.001
@GW 17 11.891 0.000
@GW 8 -14.622 0.000
@GW 9 -13.641 0.004
@GW 10 -12.648 0.005
@GW 11 -10.354 0.002
@GW 12 -8.453 0.002
@GW 13 3.960 0.003
@GW 14 3.965 0.005
@GW 15 5.443 0.002
@GW 16 10.094 0.004
@GW 17 12.237 0.004
Beta Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 7 -15.994 0.001
@GW 8 -14.090 0.000
@GW 9 -12.312 0.000
@GW 10 -11.111 0.000
@GW 11 -9.813 0.000
@GW 12 0.782 0.000
@GW 13 3.891 0.001
@GW 14 3.844 0.000
@GW 15 5.494 0.001
@GW 16 9.576 0.001
@GW 7 -17.756 0.002
@GW 8 -14.841 0.005
@GW 9 -13.401 0.002
@GW 10 -11.674 0.002
@GW 11 -10.239 0.002
@GW 12 1.142 0.002
@GW 13 4.416 0.002
@GW 14 4.181 0.005
@GW 15 5.574 0.003
@GW 16 10.343 0.017 *
* Result did not converge
GW Timing Statistics
------------------------
Init. : 0.1s
ERIs : 0.2s
Init. : 0.0s
ERIs : 0.1s
Distr.: 0.0s
Vxc : 0.7s
Vxc : 0.8s
Vmo : 0.0s
iW : 0.7s
iW : 0.5s
I_n : 0.0s
R_n : 0.3s
GW total: 2.0s
GW total: 1.8s
Task times cpu: 12.2s wall: 12.3s
NWChem Input Module
-------------------
NWChem GW@DFT Module
--------------------
GW Trifluoroacetamide
NWChem DFT Module
-----------------
GW Trifluoroacetamide
Summary of "ao basis" -> "ao basis" (spherical)
------------------------------------------------------------------------------
Tag Description Shells Functions and Types
---------------- ------------------------------ ------ ---------------------
O cc-pvdz 6 14 3s2p1d
N cc-pvdz 6 14 3s2p1d
C cc-pvdz 6 14 3s2p1d
H cc-pvdz 3 5 2s1p
Caching 1-el integrals
General Information
-------------------
SCF calculation type: DFT
Wavefunction type: spin polarized.
No. of atoms : 5
No. of electrons : 23
Alpha electrons : 12
Beta electrons : 11
Charge : 0
Spin multiplicity: 2
Use of symmetry is: off; symmetry adaption is: off
Maximum number of iterations: 50
This is a Direct SCF calculation.
AO basis - number of functions: 52
number of shells: 24
A Charge density fitting basis will be used.
CD basis - number of functions: 196
number of shells: 66
Convergence on energy requested: 1.00D-08
Convergence on density requested: 1.00D-05
Convergence on gradient requested: 5.00D-04
XC Information
--------------
PBE96 Method XC Functional
PerdewBurkeErnzerhof Exchange Functional 1.000
Perdew 1991 LDA Correlation Functional 1.000 local
PerdewBurkeErnz. Correlation Functional 1.000 non-local
Grid Information
----------------
Grid used for XC integration: fine
Radial quadrature: Mura-Knowles
Angular quadrature: Lebedev.
Tag B.-S. Rad. Rad. Pts. Rad. Cut. Ang. Pts.
--- ---------- --------- --------- ---------
O 0.60 70 10.0 590
N 0.65 70 7.0 590
C 0.70 70 8.0 590
H 0.35 60 8.0 590
Grid pruning is: on
Number of quadrature shells: 330
Spatial weights used: Erf1
Convergence Information
-----------------------
Convergence aids based upon iterative change in
total energy or number of iterations.
Levelshifting, if invoked, occurs when the
HOMO/LUMO gap drops below (HL_TOL): 1.00D-02
DIIS, if invoked, will attempt to extrapolate
using up to (NFOCK): 10 stored Fock matrices.
Damping( 0%) Levelshifting(0.5) DIIS
--------------- ------------------- ---------------
dE on: start ASAP start
dE off: 2 iters 50 iters 50 iters
Screening Tolerance Information
-------------------------------
Density screening/tol_rho: 1.00D-11
AO Gaussian exp screening on grid/accAOfunc: 18
CD Gaussian exp screening on grid/accCDfunc: 20
XC Gaussian exp screening on grid/accXCfunc: 20
Schwarz screening/accCoul: 1.00D-14
Loading old vectors from job with title :
GW Trifluoroacetamide
Time after variat. SCF: 12.5
3 Center 2 Electron Integral Information
----------------------------------------
Maximum number of 3-center 2e- integrals is: 529984.
This is reduced with Schwarz screening to: 367500.
Incore requires a per proc buffer size of: 175029.
Minimum dble words available (all nodes) is: 131065476
This is reduced (for later use) to: 130683464
proc 0 Suggested buffer size is: 175029
Max Suggested buffer size is: 175029
no. integral batches is: 1000
Incore memory use for 3-center 2e- integrals is turned off.
Time prior to 1st pass: 12.5
Total DFT energy = -169.066312547888
One electron energy = -355.393692322893
Coulomb energy = 144.424774055736
Exchange-Corr. energy = -21.405772779389
Nuclear repulsion energy = 63.308378498658
Numeric. integr. density = 22.999999954126
Total iterative time = 0.9s
alpha - beta orbital overlaps
-----------------------------
alpha 1 2 3 4 5 6 7 8 9 10
beta 1 2 3 4 5 6 7 8 9 10
overlap 1.000 1.000 1.000 1.000 1.000 0.996 0.996 0.990 1.000 0.985
alpha 11 12 13 14 15 16 17 18 19 20
beta 11 12 13 14 15 16 17 18 19 20
overlap 1.000 0.993 0.999 0.994 0.995 0.999 0.999 0.952 1.000 0.945
alpha 21 22 23 24 25 26 27 28 29 30
beta 21 22 23 24 25 26 27 28 29 30
overlap 0.958 0.999 0.998 1.000 0.998 0.999 0.998 0.996 0.997 0.980
alpha 31 32 33 34 35 36 37 38 39 40
beta 31 32 33 34 35 36 37 38 39 40
overlap 0.985 0.996 0.991 0.992 0.999 1.000 0.997 0.994 0.997 0.999
alpha 41 42 43 44 45 46 47 48 49 50
beta 41 42 43 44 45 46 47 48 49 50
overlap 0.999 0.996 0.994 0.999 0.998 0.998 1.000 0.999 0.999 1.000
alpha 51 52
beta 51 52
overlap 0.999 0.999
--------------------------
Expectation value of S2:
--------------------------
<S2> = 0.7519 (Exact = 0.7500)
center of mass
--------------
x = -0.02948626 y = 0.13290476 z = -0.00001949
moments of inertia (a.u.)
------------------
162.621484064576 0.000000000000 -0.000000000000
0.000000000000 16.578708728791 -0.000000000000
-0.000000000000 -0.000000000000 179.200070028499
Multipole analysis of the density
---------------------------------
L x y z total alpha beta nuclear
- - - - ----- ----- ---- -------
0 0 0 0 -0.000000 -12.000000 -11.000000 23.000000
1 1 0 0 0.238427 0.541429 -0.303002 0.000000
1 0 1 0 -1.321737 -0.969708 -0.352030 0.000000
1 0 0 1 0.000663 -0.000054 0.000717 0.000000
2 2 0 0 -11.884337 -13.068472 -9.699498 10.883633
2 1 1 0 -0.250985 0.506239 0.213062 -0.970286
2 1 0 1 0.006019 -0.000243 -0.000755 0.007017
2 0 2 0 -11.275818 -52.448505 -49.703057 90.875744
2 0 1 1 0.001863 0.000480 -0.000509 0.001891
2 0 0 2 -13.100006 -6.900236 -6.199822 0.000052
***********************
*** GW@DFT ***
***********************
Memory Information
------------------
Available GA per process is 1000.0 MBs
Available stack per process is 500.0 MBs
Available heap per process is 500.0 MBs
Peak memory in Initialization
Peak GA per process is 1.4 MBs
Peak stack per process is 0.0 MBs
Peak heap per process is 1.6 MBs
Memory not in MA per process 0.0 MBs
Peak memory in AFF Calculation
Peak GA per process is 5.2 MBs
Peak stack per process is 6.6 MBs
Peak heap per process is 0.0 MBs
Memory not in MA per process 6.6 MBs
Computing 2-center integrals
Computing Inverse Cholesky factor
Computing 3-center integrals
Computing V_xc ...
Computing Sigma_x
Computing RPA polarizability
Solving quasiparticle equations
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 8 -13.237 0.001
@GW 9 -12.542 0.000
@GW 10 -12.360 0.000
@GW 11 -9.991 0.000
@GW 12 -7.780 0.001
@GW 13 3.732 0.000
@GW 14 3.851 0.000
@GW 15 5.314 0.000
@GW 16 9.843 0.001
@GW 17 11.907 0.000
Beta Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 7 -16.569 0.000
@GW 8 -14.089 0.000
@GW 9 -12.853 0.000
@GW 10 -11.113 0.000
@GW 11 -9.813 0.000
@GW 12 0.782 0.000
@GW 13 4.374 0.000
@GW 14 4.112 0.000
@GW 15 5.494 0.000
@GW 16 10.421 0.000
GW Timing Statistics
------------------------
Init. : 0.1s
ERIs : 0.2s
Distr.: 0.0s
Vxc : 0.7s
Sigma_x : 0.0s
RPA : 0.3s
Wmn : 0.0s
Sigma_c : 0.0s
GW total: 1.3s
Task times cpu: 2.3s wall: 2.3s
Task times cpu: 13.8s wall: 13.8s
NWChem Input Module
@ -1160,11 +864,11 @@ MA usage statistics:
heap stack
---- -----
current number of blocks 0 0
maximum number of blocks 28 58
maximum number of blocks 28 57
current total bytes 0 0
maximum total bytes 21812064 30217960
maximum total K-bytes 21813 30218
maximum total M-bytes 22 31
maximum total bytes 23306336 30269032
maximum total K-bytes 23307 30270
maximum total M-bytes 24 31
CITATION
@ -1200,24 +904,25 @@ MA usage statistics:
AUTHORS
-------
E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, M. Valiev, W. A. de Jong,
T. P. Straatsma, H. J. J. van Dam, D. Wang, T. L. Windus, N. P. Bauman,
A. Panyala, J. Hammond, J. Autschbach, K. Bhaskaran-Nair, J. Brabec,
K. Lopata, S. A. Fischer, S. Krishnamoorthy, M. Jacquelin, W. Ma, M. Klemm,
O. Villa, Y. Chen, V. Anisimov, F. Aquino, S. Hirata, M. T. Hackler,
Eric Hermes, L. Jensen, J. E. Moore, J. C. Becca, V. Konjkov,
D. Mejia-Rodriguez, T. Risthaus, M. Malagoli, A. Marenich,
A. Otero-de-la-Roza, J. Mullin, P. Nichols, R. Peverati, J. Pittner, Y. Zhao,
P.-D. Fan, A. Fonari, M. J. Williamson, R. J. Harrison, J. R. Rehr,
M. Dupuis, D. Silverstein, D. M. A. Smith, J. Nieplocha, V. Tipparaju,
M. Krishnan, B. E. Van Kuiken, A. Vazquez-Mayagoitia, M. Swart, Q. Wu,
T. Van Voorhis, A. A. Auer, M. Nooijen, L. D. Crosby, E. Brown, G. Cisneros,
G. I. Fann, H. Fruchtl, J. Garza, K. Hirao, R. A. Kendall, J. A. Nichols,
K. Tsemekhman, K. Wolinski, J. Anchell, D. E. Bernholdt, P. Borowski,
T. Clark, D. Clerc, H. Dachsel, M. J. O. Deegan, K. Dyall, D. Elwood,
E. Glendening, M. Gutowski, A. C. Hess, J. Jaffe, B. G. Johnson, J. Ju,
R. Kobayashi, R. Kutteh, Z. Lin, R. Littlefield, X. Long, B. Meng,
T. Nakajima, S. Niu, L. Pollack, M. Rosing, K. Glaesemann, G. Sandrone,
M. Stave, H. Taylor, G. Thomas, J. H. van Lenthe, A. T. Wong, Z. Zhang.
E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, D. Mejia-Rodriguez,
N. P. Bauman, A. Panyala, R. J. Harrison, M. Valiev, W. A. de Jong,
T. P. Straatsma, H. J. J. van Dam, D. Wang, T. L. Windus, J. Hammond,
J. Autschbach, R. de P. Soares, A. Kunitsa, K. Bhaskaran-Nair, J. Brabec,
K. Lopata, S. A. Fischer, S. Krishnamoorthy, M. Jacquelin, A. Woods, W. Ma,
M. Klemm, O. Villa, Y. Chen, V. Anisimov, F. Aquino, S. Hirata,
M. T. Hackler, E. Hermes, L. Jensen, J. E. Moore, J. C. Becca, V. Konjkov,
T. Risthaus, M. Malagoli, A. Marenich, A. Otero-de-la-Roza, J. Mullin,
P. Nichols, R. Peverati, J. Pittner, Y. Zhao, P.-D. Fan, A. Fonari,
M. J. Williamson, J. R. Rehr, M. Dupuis, D. Silverstein, D. M. A. Smith,
J. Nieplocha, V. Tipparaju, M. Krishnan, B. E. Van Kuiken,
A. Vazquez-Mayagoitia, M. Swart, Q. Wu, T. Van Voorhis, A. A. Auer,
M. Nooijen, L. D. Crosby, E. Brown, G. Cisneros, G. I. Fann, H. Fruchtl,
J. Garza, K. Hirao, R. A. Kendall, J. A. Nichols, K. Tsemekhman, K. Wolinski,
J. Anchell, D. E. Bernholdt, P. Borowski, T. Clark, D. Clerc, H. Dachsel,
M. J. O. Deegan, K. Dyall, D. Elwood, E. Glendening, M. Gutowski, A. C. Hess,
J. Jaffe, B. G. Johnson, J. Ju, R. Kobayashi, R. Kutteh, Z. Lin,
R. Littlefield, X. Long, B. Meng, T. Nakajima, S. Niu, L. Pollack, M. Rosing,
K. Glaesemann, G. Sandrone, M. Stave, H. Taylor, G. Thomas, J. H. van Lenthe,
A. T. Wong, Z. Zhang.
Total times cpu: 14.5s wall: 14.7s
Total times cpu: 14.2s wall: 14.2s

View file

@ -1,5 +1,5 @@
argument 1 = gw_symmetry.nw
NWChem w/ OpenMP: maximum threads = 8
argument 1 = /people/meji656/calculations/bse/QA/gw_symmetry.nw
NWChem w/ OpenMP: maximum threads = 1
@ -53,7 +53,7 @@ task dft gw
Northwest Computational Chemistry Package (NWChem) 7.2.0
Northwest Computational Chemistry Package (NWChem) 7.3.0
--------------------------------------------------------
@ -61,7 +61,7 @@ task dft gw
Pacific Northwest National Laboratory
Richland, WA 99352
Copyright (c) 1994-2022
Copyright (c) 1994-2025
Pacific Northwest National Laboratory
Battelle Memorial Institute
@ -86,21 +86,21 @@ task dft gw
Job information
---------------
hostname = WE41476
program = nwchem
date = Wed Jul 12 13:50:26 2023
hostname = dc005
program = /qfs/people/meji656/software/nwchem.bse/bin/LINUX64/nwchem
date = Thu Oct 16 10:35:32 2025
compiled = Wed_Jul_12_13:50:03_2023
source = /Users/meji656/Sources/nwchem_bse
nwchem branch = 7.2.0
nwchem revision = nwchem_on_git-4850-g011c901253
ga revision = 5.8.0
use scalapack = T
input = gw_symmetry.nw
compiled = Thu_Oct_16_10:07:15_2025
source = /qfs/people/meji656/software/nwchem.bse
nwchem branch = 7.3.0
nwchem revision = v7.2.0-beta1-1624-g0dc011d0a0
ga revision = 5.9.2
use scalapack = F
input = /people/meji656/calculations/bse/QA/gw_symmetry.nw
prefix = gw_symmetry.
data base = ./gw_symmetry.db
status = startup
nproc = 1
nproc = 7
time left = -1s
@ -108,10 +108,10 @@ task dft gw
Memory information
------------------
heap = 65535996 doubles = 500.0 Mbytes
stack = 65536001 doubles = 500.0 Mbytes
global = 131072000 doubles = 1000.0 Mbytes (distinct from heap & stack)
total = 262143997 doubles = 2000.0 Mbytes
heap = 65535998 doubles = 500.0 Mbytes
stack = 65535995 doubles = 500.0 Mbytes
global = 131072000 doubles = 1000.0 Mbytes (distinct from heap & stack)
total = 262143993 doubles = 2000.0 Mbytes
verify = yes
hardfail = no
@ -462,33 +462,32 @@ task dft gw
6 a1 7 t2 8 t2 9 t2 10 t2
11 t2 12 t2 13 a1 14 t2 15 t2
WARNING: movecs_in_org=atomic not equal to movecs_in=./gw_symmetry.movecs
Time after variat. SCF: 0.3
Time prior to 1st pass: 0.3
Time after variat. SCF: 0.4
Time prior to 1st pass: 0.4
Memory utilization after 1st SCF pass:
Heap Space remaining (MW): 65.53 65532636
Stack Space remaining (MW): 65.54 65535396
Heap Space remaining (MW): 65.53 65532638
Stack Space remaining (MW): 65.54 65535388
convergence iter energy DeltaE RMS-Dens Diis-err time
---------------- ----- ----------------- --------- --------- --------- ------
d= 0,ls=0.0,diis 1 -40.4526766257 -5.40D+01 2.73D-03 1.58D-01 0.6
d= 0,ls=0.0,diis 2 -40.4605212533 -7.84D-03 1.18D-03 2.13D-02 0.8
d= 0,ls=0.0,diis 3 -40.4604458885 7.54D-05 6.05D-04 2.16D-02 1.1
d= 0,ls=0.0,diis 4 -40.4630635417 -2.62D-03 2.24D-05 6.66D-06 1.3
d= 0,ls=0.0,diis 5 -40.4630646869 -1.15D-06 2.13D-06 1.39D-07 1.6
d= 0,ls=0.0,diis 6 -40.4630647071 -2.03D-08 7.01D-08 1.45D-10 1.9
d= 0,ls=0.0,diis 1 -40.4526766257 -5.40D+01 2.73D-03 1.58D-01 0.5
d= 0,ls=0.0,diis 2 -40.4605212533 -7.84D-03 1.18D-03 2.13D-02 0.5
d= 0,ls=0.0,diis 3 -40.4604458885 7.54D-05 6.05D-04 2.16D-02 0.6
d= 0,ls=0.0,diis 4 -40.4630635417 -2.62D-03 2.24D-05 6.66D-06 0.6
d= 0,ls=0.0,diis 5 -40.4630646869 -1.15D-06 2.13D-06 1.39D-07 0.7
d= 0,ls=0.0,diis 6 -40.4630647071 -2.03D-08 7.01D-08 1.45D-10 0.7
Total DFT energy = -40.463064707119
One electron energy = -79.987919842739
Coulomb energy = 32.854590174254
Total DFT energy = -40.463064707118
One electron energy = -79.987919842738
Coulomb energy = 32.854590174255
Exchange-Corr. energy = -6.837400210970
Nuclear repulsion energy = 13.507665172336
Numeric. integr. density = 10.000008758007
Total iterative time = 1.5s
Total iterative time = 0.3s
@ -523,11 +522,11 @@ task dft gw
1 1 0 0 -0.000000 -0.000000 -0.000000 0.000000
1 0 1 0 0.000000 0.000000 0.000000 0.000000
1 0 0 1 -0.000000 -0.000000 -0.000000 0.000000
1 0 0 1 0.000000 0.000000 0.000000 0.000000
2 2 0 0 -6.283208 -5.939937 -5.939937 5.596666
2 1 1 0 0.000000 0.000000 0.000000 0.000000
2 1 0 1 0.000000 0.000000 0.000000 0.000000
2 1 1 0 -0.000000 -0.000000 -0.000000 0.000000
2 1 0 1 -0.000000 -0.000000 -0.000000 0.000000
2 0 2 0 -6.283208 -5.939937 -5.939937 5.596666
2 0 1 1 -0.000000 -0.000000 -0.000000 0.000000
2 0 0 2 -6.283208 -5.939937 -5.939937 5.596666
@ -626,18 +625,18 @@ task dft gw
Memory Information
------------------
Available GA per process is 999.9 MBs
Available GA per process is 1000.0 MBs
Available stack per process is 500.0 MBs
Available heap per process is 499.9 MBs
Peak memory in Initialization
Peak GA per process is 1.1 MBs
Peak GA per process is 0.2 MBs
Peak stack per process is 0.1 MBs
Peak heap per process is 1.2 MBs
Peak heap per process is 0.5 MBs
Memory not in MA per process 0.0 MBs
Peak memory in AFF Calculation
Peak GA per process is 2.5 MBs
Peak GA per process is 0.4 MBs
Peak stack per process is 1.3 MBs
Peak heap per process is 0.0 MBs
Memory not in MA per process 1.3 MBs
@ -657,25 +656,25 @@ task dft gw
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 5 -13.866 0.000
@GW 5 -13.864 0.003
Applying -4.4168 eV scissor-shift to rest of particle states
Applying -4.4144 eV scissor-shift to rest of particle states
GW Timing Statistics
------------------------
Init. : 0.1s
ERIs : 0.3s
ERIs : 0.1s
Distr.: 0.0s
Vxc : 0.2s
Vxc : 0.0s
Sigma_x : 0.0s
RPA : 0.0s
Wmn : 0.0s
Sigma_c : 0.0s
GW total: 0.5s
GW total: 0.2s
Task times cpu: 2.4s wall: 2.4s
Task times cpu: 0.6s wall: 0.6s
NWChem Input Module
@ -784,7 +783,7 @@ task dft gw
Loading old vectors from job with title :
GW Methane Td
W vals
Symmetry analysis of molecular orbitals - initial
@ -800,28 +799,28 @@ GW Methane Td
6 a1 7 t2 8 t2 9 t2 10 t2
11 t2 12 t2 13 t2 14 t2 15 t2
Time after variat. SCF: 2.5
Time prior to 1st pass: 2.5
Time after variat. SCF: 0.9
Time prior to 1st pass: 0.9
Memory utilization after 1st SCF pass:
Heap Space remaining (MW): 65.53 65532636
Stack Space remaining (MW): 65.54 65535396
Heap Space remaining (MW): 65.53 65532638
Stack Space remaining (MW): 65.54 65535388
convergence iter energy DeltaE RMS-Dens Diis-err time
---------------- ----- ----------------- --------- --------- --------- ------
d= 0,ls=0.0,diis 1 -40.4630647077 -5.40D+01 1.34D-08 1.68D-12 3.5
d= 0,ls=0.0,diis 2 -40.4630647077 -2.15D-12 7.98D-09 3.73D-12 4.4
d= 0,ls=0.0,diis 1 -40.4630647076 -5.40D+01 1.34D-08 1.68D-12 1.0
d= 0,ls=0.0,diis 2 -40.4630647076 2.13D-13 7.98D-09 3.73D-12 1.1
Total DFT energy = -40.463064707739
One electron energy = -79.987938823681
Coulomb energy = 32.854611046058
Exchange-Corr. energy = -6.837402102452
Total DFT energy = -40.463064707584
One electron energy = -79.987938823837
Coulomb energy = 32.854611046046
Exchange-Corr. energy = -6.837402102129
Nuclear repulsion energy = 13.507665172336
Numeric. integr. density = 10.000008760592
Numeric. integr. density = 10.000008760422
Total iterative time = 1.9s
Total iterative time = 0.3s
@ -854,15 +853,15 @@ GW Methane Td
- - - - ----- ----- ---- -------
0 0 0 0 -0.000000 -5.000000 -5.000000 10.000000
1 1 0 0 -0.000000 -0.000000 -0.000000 0.000000
1 0 1 0 -0.000000 -0.000000 -0.000000 0.000000
1 0 0 1 -0.000000 -0.000000 -0.000000 0.000000
1 1 0 0 0.000000 0.000000 0.000000 0.000000
1 0 1 0 0.000000 0.000000 0.000000 0.000000
1 0 0 1 0.000000 0.000000 0.000000 0.000000
2 2 0 0 -6.283210 -5.939938 -5.939938 5.596666
2 1 1 0 -0.000000 -0.000000 -0.000000 0.000000
2 1 0 1 -0.000000 -0.000000 -0.000000 0.000000
2 1 1 0 0.000000 0.000000 0.000000 0.000000
2 1 0 1 0.000000 0.000000 0.000000 0.000000
2 0 2 0 -6.283210 -5.939938 -5.939938 5.596666
2 0 1 1 -0.000000 -0.000000 -0.000000 0.000000
2 0 1 1 0.000000 0.000000 0.000000 0.000000
2 0 0 2 -6.283210 -5.939938 -5.939938 5.596666
@ -873,18 +872,18 @@ GW Methane Td
Memory Information
------------------
Available GA per process is 999.9 MBs
Available GA per process is 1000.0 MBs
Available stack per process is 500.0 MBs
Available heap per process is 499.9 MBs
Peak memory in Initialization
Peak GA per process is 1.1 MBs
Peak GA per process is 0.2 MBs
Peak stack per process is 0.1 MBs
Peak heap per process is 1.2 MBs
Peak heap per process is 0.5 MBs
Memory not in MA per process 0.0 MBs
Peak memory in AFF Calculation
Peak GA per process is 2.5 MBs
Peak GA per process is 0.4 MBs
Peak stack per process is 1.3 MBs
Peak heap per process is 0.0 MBs
Memory not in MA per process 1.3 MBs
@ -904,25 +903,25 @@ GW Methane Td
Alpha Orbitals
State Energy (eV) Error (eV)
-----------------------------------
@GW 5 -13.866 0.000
@GW 5 -13.864 0.003
Applying -4.4168 eV scissor-shift to rest of particle states
Applying -4.4144 eV scissor-shift to rest of particle states
GW Timing Statistics
------------------------
Init. : 0.0s
ERIs : 0.2s
ERIs : 0.0s
Distr.: 0.0s
Vxc : 0.5s
Vxc : 0.1s
Sigma_x : 0.0s
RPA : 0.0s
Wmn : 0.0s
Sigma_c : 0.0s
GW total: 0.8s
GW total: 0.2s
Task times cpu: 2.8s wall: 2.8s
Task times cpu: 0.5s wall: 0.5s
NWChem Input Module
@ -944,9 +943,9 @@ MA usage statistics:
current number of blocks 0 0
maximum number of blocks 29 62
current total bytes 0 0
maximum total bytes 32636240 32078504
maximum total K-bytes 32637 32079
maximum total M-bytes 33 33
maximum total bytes 31929024 32078520
maximum total K-bytes 31930 32079
maximum total M-bytes 32 33
CITATION
@ -982,24 +981,25 @@ MA usage statistics:
AUTHORS
-------
E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, M. Valiev, D. Mejia-Rodriguez,
A. Kunitsa, N. P. Bauman, A. Panyala, W. A. de Jong, T. P. Straatsma,
H. J. J. van Dam, D. Wang, T. L. Windus, J. Hammond, J. Autschbach, A. Woods,
K. Bhaskaran-Nair, J. Brabec, K. Lopata, S. A. Fischer, S. Krishnamoorthy,
M. Jacquelin, W. Ma, M. Klemm, O. Villa, Y. Chen, V. Anisimov, F. Aquino,
E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, D. Mejia-Rodriguez,
N. P. Bauman, A. Panyala, R. J. Harrison, M. Valiev, W. A. de Jong,
T. P. Straatsma, H. J. J. van Dam, D. Wang, T. L. Windus, J. Hammond,
J. Autschbach, R. de P. Soares, A. Kunitsa, K. Bhaskaran-Nair, J. Brabec,
K. Lopata, S. A. Fischer, S. Krishnamoorthy, M. Jacquelin, R. Goswami,
A. Woods, W. Ma, M. Klemm, O. Villa, Y. Chen, V. Anisimov, F. Aquino,
S. Hirata, M. T. Hackler, E. Hermes, L. Jensen, J. E. Moore, J. C. Becca,
V. Konjkov, T. Risthaus, M. Malagoli, A. Marenich, A. Otero-de-la-Roza,
J. Mullin, P. Nichols, R. Peverati, J. Pittner, Y. Zhao, P.-D. Fan,
A. Fonari, M. J. Williamson, R. J. Harrison, J. R. Rehr, M. Dupuis,
D. Silverstein, D. M. A. Smith, J. Nieplocha, V. Tipparaju, M. Krishnan,
B. E. Van Kuiken, A. Vazquez-Mayagoitia, M. Swart, Q. Wu, T. Van Voorhis,
A. A. Auer, M. Nooijen, L. D. Crosby, E. Brown, G. Cisneros, G. I. Fann,
H. Fruchtl, J. Garza, K. Hirao, R. A. Kendall, J. A. Nichols, K. Tsemekhman,
K. Wolinski, J. Anchell, D. E. Bernholdt, P. Borowski, T. Clark, D. Clerc,
H. Dachsel, M. J. O. Deegan, K. Dyall, D. Elwood, E. Glendening, M. Gutowski,
A. C. Hess, J. Jaffe, B. G. Johnson, J. Ju, R. Kobayashi, R. Kutteh, Z. Lin,
A. Fonari, M. J. Williamson, J. R. Rehr, M. Dupuis, D. Silverstein,
D. M. A. Smith, J. Nieplocha, V. Tipparaju, M. Krishnan, B. E. Van Kuiken,
A. Vazquez-Mayagoitia, M. Swart, Q. Wu, T. Van Voorhis, A. A. Auer,
M. Nooijen, L. D. Crosby, E. Brown, G. Cisneros, G. I. Fann, H. Fruchtl,
J. Garza, K. Hirao, R. A. Kendall, J. A. Nichols, K. Tsemekhman, K. Wolinski,
J. Anchell, D. E. Bernholdt, P. Borowski, T. Clark, D. Clerc, H. Dachsel,
M. J. O. Deegan, K. Dyall, D. Elwood, E. Glendening, M. Gutowski, A. C. Hess,
J. Jaffe, B. G. Johnson, J. Ju, R. Kobayashi, R. Kutteh, Z. Lin,
R. Littlefield, X. Long, B. Meng, T. Nakajima, S. Niu, L. Pollack, M. Rosing,
K. Glaesemann, G. Sandrone, M. Stave, H. Taylor, G. Thomas, J. H. van Lenthe,
A. T. Wong, Z. Zhang.
Total times cpu: 5.2s wall: 5.2s
Total times cpu: 1.3s wall: 1.3s

44
src/bsemol/GNUmakefile Normal file
View file

@ -0,0 +1,44 @@
OBJ = bsemol.o bse_input.o bse_init.o bse_analytic.o \
bse_davidson.o bse_lanczos.o bse_defaults.o \
bse_finalize.o bse_ri_init.o bse_wmn.o \
bse_analytic_cis.o bse_analytic_tda.o \
bse_oscstr.o bse_puteris.o bse_get_eia.o \
bse_davidson_diagonal.o \
bse_davidson_kvec.o bse_davidson_guess.o \
bse_davidson_kortho.o bse_davidson_restart.o \
bse_davidson_newvecs.o bse_davidson_residual.o \
bse_s2.o
USES_BLAS = bse_buildw.F bse_ri_init.F bse_wmn.F \
bse_davidson.F bse_davidson_diagonal.F bse_davidson_kvec.F bse_davidson_kortho.F \
bse_init.F bse_puteris.F bse_s2.F
OBJ_OPTIMIZE = bse_buildw.o
LIBRARY = libbsemol.a
HEADERS = ../nwdft/include/cdft.fh ../nwdft/grid/cgridfile.fh ../nwdft/include/dft_scaleig.fh
LIB_DEFINES =
LIB_INCLUDES = -I./
include ../config/makefile.h
ifeq ($(_USE_SCALAPACK),Y)
_USE_ELPA = $(shell ${GA_PATH}/bin/ga-config --use_elpa| awk ' /1/ {print "Y"}')
ifeq ($(_USE_ELPA),Y)
LIB_DEFINES += -DUSE_ELPA
endif
#assuming we have GAPZHEEVD
GOTGAPZHEEVD= Y
ifeq ($(GOTGAPZHEEVD),Y)
LIB_DEFINES += -DUSE_GAPZHEEVD
endif
_PDSYEVR = $(shell [ -e ${NWCHEM_TOP}/src/tools/install/lib/libga.a ] && nm ${NWCHEM_TOP}/src/tools/install/lib/libga.a|egrep T |egrep ga_pdsyevr |awk '/ga_pdsyevr/ {print "Y"}')
ifeq ($(_PDSYEVR),Y)
LIB_DEFINES += -DUSE_PDSYEVR
endif
endif
include ../config/makelib.h

44
src/bsemol/bse.fh Normal file
View file

@ -0,0 +1,44 @@
double precision, parameter :: ha2ev = 27.211386245988d0
type :: bse_params_t
sequence
integer :: rtdb
integer :: me, nprocs, iMaxThreads
integer :: ipol
logical :: main, odft
c
integer :: maxiter, nspace
integer :: natoms
integer :: nmo, nri
integer :: nocc(2),nvir(2),npoles(2),mynpoles(2)
integer :: oolo(2),oohi(2)
integer :: ovlo(2),ovhi(2)
integer :: vvlo(2),vvhi(2)
integer :: llmo(2),nqp(2)
integer :: nroots,nmax
double precision :: tol2e, eta
c
integer :: l_mf_evals, k_mf_evals
integer :: l_gw_evals, k_gw_evals
integer :: l_occ, k_occ
integer :: l_irs(2),k_irs(2)
integer :: l_map,k_map
integer :: nblocks
c
integer :: fd
c
integer :: g_movecs(2), g_moeris, g_sols
integer :: g_erioo(2), g_eriov(2), g_erivv(2), g_erim
integer :: g_wov(2), g_apb, g_dipa(3), g_dipb(3), g_amb
integer :: g_quada(6), g_quadb(6)
integer :: g_dipmaga(3), g_dipmagb(3)
logical :: noio
c
logical :: ncap, davidson, lanczos, analytic, tda
logical :: singlet, triplet, dowindow
double precision :: thresh, fdstep, vxddp(2), vxddm(2)
double precision :: elower, eupper
end type bse_params_t

82
src/bsemol/bse_analytic.F Normal file
View file

@ -0,0 +1,82 @@
subroutine bse_analytic(pars)
implicit none
#include "global.fh"
#include "errquit.fh"
#include "mafdecls.fh"
#include "bse.fh"
#include "cdft.fh"
type(bse_params_t) :: pars
character(*), parameter :: pname = 'bse_analytic: '
integer lSize,isp
integer l_omega, k_omega
integer l_dia,k_dia
integer lW,kW
integer maxpoles,maxnpoles,totnpoles,nri
! synchronize parallel processes
call ga_sync
!initialize useful variables
nri = pars%nri
maxpoles = maxval(pars%mynpoles(1:pars%ipol))
maxnpoles = maxval(pars%npoles(1:pars%ipol))
totnpoles = maxnpoles*pars%ipol
! allocation
if(.not.ma_push_get(mt_dbl,totnpoles,'omega',l_omega,k_omega))
& call errquit(pname//'failed to allocate omega',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,totnpoles,'dia',l_dia,k_dia))
& call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,nri**2,'W',lW,kW))
& call errquit(pname//'failed to allocate W',0,MA_ERR)
! obtain "old" eigenvalue difference to build W
do isp=1,pars%ipol
lSize = (isp-1)*maxpoles
call bse_get_eia(dbl_mb(pars%k_mf_evals+(isp-1)*nbf_ao),
& dbl_mb(k_dia+lSize),pars%nocc(isp),
& pars%nvir(isp),pars%ovlo(isp),pars%ovhi(isp))
enddo
! obtain screened Coulomb matrix
call bse_buildw(pars,dbl_mb(k_dia),dbl_mb(kW),pars%nmo,nri,
& maxpoles,pars%ipol)
! transform ERIs using screened Coulomb matrix
call bse_wmn(pars,dbl_mb(kW),pars%nmo,nri,pars%ipol)
if(.not.ma_chop_stack(lW))
& call errquit(pname//'failed to chop stack',0,MA_ERR)
! obtain "GW" eigenvalue difference to build Casida matrix
do isp=1,pars%ipol
lSize = (isp-1)*maxpoles
call bse_get_eia(dbl_mb(pars%k_gw_evals+(isp-1)*nbf_ao),
& dbl_mb(k_dia+lSize),pars%nocc(isp),
& pars%nvir(isp),pars%ovlo(isp),
& pars%ovhi(isp))
enddo
! Build and diagonalize Casida-like matrix
if (pars%tda) then
call bse_analytic_tda(pars,dbl_mb(k_dia),dbl_mb(k_omega),
& maxpoles,totnpoles,pars%singlet)
else
call bse_analytic_cis(pars,dbl_mb(k_dia),dbl_mb(k_omega),
& maxpoles,totnpoles,pars%singlet)
endif
! Compute oscillator strengths
call bse_oscstr(pars,dbl_mb(k_omega),totnpoles)
! Destroy GAs
if (.not.ga_destroy(pars%g_apb))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
! chop stack
if(.not.ma_chop_stack(l_omega))
$ call errquit(pname//'failed to chop stack',0,MA_ERR)
end subroutine bse_analytic

View file

@ -0,0 +1,359 @@
subroutine bse_analytic_cis(pars,wia,omega,maxpoles,totpoles,
$ singlet)
implicit none
#include "mafdecls.fh"
#include "errquit.fh"
#include "global.fh"
#include "bse.fh"
#ifdef SCALAPACK
#include "scaleig.fh"
integer ga_cholesky,ga_llt_i
external ga_cholesky,ga_llt_i
#endif
type(bse_params_t) :: pars
logical singlet
integer maxpoles, totpoles, npoles
double precision wia(maxpoles,pars%ipol)
double precision omega(totpoles)
character(*),parameter :: pname = 'bse_casida: '
integer isp,jsp,nri,off1,off2
integer ilo,ihi,jlo,jhi
integer nocc2(2),nvir2(2)
integer ipole,imo,jmo,amo,bmo
integer k,l,kproc,lproc,klocal,llocal,kglobal,lglobal
integer g_diag, g_tmp, g_amb, ktmp, ld, info, i, j
integer ij,jb,llia,ulia
#ifdef USE_OPENMP
call util_blas_set_num_threads(pars%iMaxthreads)
#endif
nri = pars%nri
npoles = pars%npoles(1)
if (pars%ipol.gt.1) npoles = npoles + pars%npoles(2)
if(.not.ga_create(mt_dbl,npoles,npoles,'apb',0,npoles,
& pars%g_apb))
& call errquit(pname//'could not create A+B GA',0,GA_ERR)
if(.not.ga_create(mt_dbl,npoles,npoles,'amb',0,npoles,
& g_amb))
& call errquit(pname//'could not create A-B GA',0,GA_ERR)
if (.not.nga_create(mt_dbl,1,npoles,'wia',0,g_diag))
$ call errquit(pname//'could not create wia',0,GA_ERR)
if (pars%ipol.eq.1) then
nocc2(1) = pars%nocc(1)**2
nvir2(1) = pars%nvir(1)**2
call ga_distribution(pars%g_apb,pars%me,ilo,ihi,jlo,jhi)
if (jhi.ge.jlo) call nga_put(g_diag,ilo,ihi,wia,ihi-ilo+1)
! Screened Coulomb contribution from A
if(.not.ga_create(mt_dbl,nvir2(1),nocc2(1),'tmp',nvir2(1),0,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('t','n',nvir2(1),nocc2(1),pars%nri,-1d0,
% pars%g_erivv(1),pars%g_erioo(1),0d0,g_tmp)
call ga_distribution(g_tmp,pars%me,ilo,ihi,jlo,jhi)
call ga_access(g_tmp,ilo,ihi,jlo,jhi,ktmp,ld)
! Need to take care about indexing
do ij=jlo,jhi
jmo = (ij-1)/pars%nocc(1) + 1
imo = ij - (jmo-1)*pars%nocc(1)
llia = (imo-1)*pars%nvir(1) + 1
ulia = imo*pars%nvir(1)
do bmo = 1,pars%nvir(1)
jb = (jmo-1)*pars%nvir(1) + bmo
call ga_put(pars%g_apb,llia,ulia,jb,jb,
& dbl_mb(ktmp+(ij-jlo)*nvir2(1)+(bmo-1)*pars%nvir(1)),
& pars%nvir(1))
enddo
enddo
call ga_sync()
call ga_copy(pars%g_apb,g_amb)
call ga_release(g_tmp,ilo,ihi,jlo,jhi)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
! Screened-Coulomb contribution from B
if(.not.ga_create(mt_dbl,npoles,npoles,'tmp',0,npoles,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('t','n',npoles,npoles,pars%nri,-1d0,
& pars%g_wov(1),pars%g_wov(1),0d0,g_tmp)
call ga_add(1d0,g_tmp,1d0,pars%g_apb,pars%g_apb)
call ga_add(-1d0,g_tmp,1d0,g_amb,g_amb)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
if (singlet) then
! Hartree contribution
call ga_dgemm('t','n',npoles,npoles,pars%nri,4d0,
& pars%g_eriov(1),pars%g_eriov(1),1d0,pars%g_apb)
endif
! Diagonal term
call ga_add_diagonal(pars%g_apb,g_diag)
call ga_add_diagonal(g_amb,g_diag)
! Destroy ERIs
if (.not.ga_destroy(pars%g_erioo(1)))
& call errquit(pname//'could not destroy erioo GA',1,GA_ERR)
if (.not.ga_destroy(pars%g_erivv(1)))
& call errquit(pname//'could not destroy erivv GA',1,GA_ERR)
if (.not.ga_destroy(pars%g_eriov(1)))
& call errquit(pname//'could not destroy eriov GA',1,GA_ERR)
if (.not.ga_destroy(pars%g_wov(1)))
& call errquit(pname//'could not destroy wov GA',1,GA_ERR)
! First step is to get (A-B)^(1/2)
call ga_sync()
#ifdef SCALAPACK
info= ga_cholesky('L',g_amb)
#elif defined(PEIGS)
call ga_chol(g_amb, g_amb, info)
#else
call ga_chol_seq(g_amb, g_amb, info)
#endif
! Transform (A+B) matrix with Cholesky factors
if(.not.ga_create(mt_dbl,npoles,npoles,'tmp',0,npoles,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('n','n',npoles,npoles,npoles,1d0,
& pars%g_apb,g_amb,0d0,g_tmp)
call ga_dgemm('t','n',npoles,npoles,npoles,1d0,
& g_amb,g_tmp,0d0,pars%g_apb)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
! Diagonalize (A+B) matrix
#ifdef PARALLEL_DIAG
#ifdef SCALAPACK
call dft_scaleig(pars%g_apb,pars%g_apb,omega,0)
#else
call ga_diag_std(pars%g_apb,pars%g_apb,omega)
#endif
#else
call ga_diag_std_seq(pars%g_apb,pars%g_apb,omega)
#endif
!
! Since we transformed the eigenvalue equation to an Hermitian
! form, we obtain the square of the neutral excitations
! (Equation (32)).
!
!$omp parallel do
do ipole=1,npoles
omega(ipole) = dsqrt(omega(ipole))
enddo
!$omp end parallel do
! Transform eigenvectors with (A-B)^(1/2)
if(.not.ga_create(mt_dbl,npoles,npoles,'tmp',0,npoles,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('n','n',npoles,npoles,npoles,1d0,
& g_amb,pars%g_apb,0d0,g_tmp)
! Gather neutral excitations vectors
call nga_distribution(g_diag,pars%me,ilo,ihi)
if(ihi.ge.ilo) then
call nga_put(g_diag,ilo,ihi,1d0/dsqrt(omega(ilo:ihi)),
& ihi-ilo+1)
endif
call ga_sync()
call ga_scale_cols(g_tmp,g_diag)
! Destroy unused GAs
call ga_copy(g_tmp,pars%g_apb)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
if (.not.ga_destroy(g_amb))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
if (.not.ga_destroy(g_diag))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
else
nocc2(1) = pars%nocc(1)**2
nvir2(1) = pars%nvir(1)**2
nocc2(2) = pars%nocc(2)**2
nvir2(2) = pars%nvir(2)**2
call ga_zero(pars%g_apb)
! Gather eigenvalue differences
call ga_distribution(pars%g_eriov(1),pars%me,ilo,ihi,jlo,jhi)
call nga_put(g_diag,jlo,jhi,wia(1,1),jhi-jlo+1)
call ga_distribution(pars%g_eriov(2),pars%me,ilo,ihi,jlo,jhi)
call nga_put(g_diag,jlo+pars%npoles(1),jhi+pars%npoles(1),
$ wia(1,2),jhi-jlo+1)
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
! Add Screened Coulomb contribution from A
if (.not.ga_create(mt_dbl,nvir2(isp),nocc2(isp),'tmp',
$ nvir2(isp),0,g_tmp))
$ call errquit('bse_analytic_tda: could not create tmp',145,
$ GA_ERR)
call ga_dgemm('t','n',nvir2(isp),nocc2(isp),nri,-1d0,
$ pars%g_erivv(isp),pars%g_erioo(isp),0d0,g_tmp)
call ga_distribution(g_tmp,pars%me,ilo,ihi,jlo,jhi)
call ga_access(g_tmp,ilo,ihi,jlo,jhi,ktmp,ihi-ilo+1)
do ij=jlo,jhi
jmo = (ij-1)/pars%nocc(isp) + 1
imo = ij - (jmo-1)*pars%nocc(isp)
do bmo = 1,pars%nvir(isp)
llia = (imo-1)*pars%nvir(isp) + 1 + off1
ulia = imo*pars%nvir(isp) + off1
jb = (jmo-1)*pars%nvir(isp) + bmo + off1
call ga_put(pars%g_apb,llia,ulia,jb,jb,
& dbl_mb(ktmp+(ij-jlo)*nvir2(isp)+
& (bmo-1)*pars%nvir(isp)),pars%nvir(isp))
enddo
enddo
call ga_release(g_tmp,ilo,ihi,jlo,jhi)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',172,GA_ERR)
! Destroy ERIs
if (.not.ga_destroy(pars%g_erioo(isp)))
& call errquit(pname//'could not destroy erioo GA',175,GA_ERR)
if (.not.ga_destroy(pars%g_erivv(isp)))
& call errquit(pname//'could not destroy erivv GA',177,GA_ERR)
call ga_copy_patch('n',pars%g_apb,1+off1,
& pars%npoles(isp)+off1,1+off1,pars%npoles(isp)+off1,
& g_amb,1+off1,pars%npoles(isp)+off1,1+off1,
& pars%npoles(isp)+off1)
! Add Screened Coulomb contribution from B
if(.not.ga_create(mt_dbl,pars%npoles(isp),pars%npoles(isp),
& 'tmp',0,pars%npoles(isp),g_tmp))
& call errquit(pname//'could not create temp GA',234,GA_ERR)
call ga_dgemm('t','n',pars%npoles(isp),pars%npoles(isp),
& pars%nri,-1d0,pars%g_wov(isp),pars%g_wov(isp),
& 0d0,g_tmp)
call ga_add_patch(1d0,g_tmp,1,pars%npoles(isp),
& 1,pars%npoles(isp),
& 1d0,pars%g_apb,1+off1,pars%npoles(isp)+off1,
& 1+off1,pars%npoles(isp)+off1,
& pars%g_apb,1+off1,pars%npoles(isp)+off1,
& 1+off1,pars%npoles(isp)+off1)
call ga_add_patch(-1d0,g_tmp,1,pars%npoles(isp),
& 1,pars%npoles(isp),
& 1d0,g_amb,1+off1,pars%npoles(isp)+off1,
& 1+off1,pars%npoles(isp)+off1,
& g_amb,1+off1,pars%npoles(isp)+off1,
& 1+off1,pars%npoles(isp)+off1)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
if (.not.ga_destroy(pars%g_wov(isp)))
& call errquit(pname//'could not destroy wov GA',1,GA_ERR)
enddo
! Form RPA matrix
if (singlet) then
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
do jsp=1,pars%ipol
off2 = (jsp-1)*pars%npoles(1)
call ga_matmul_patch('t','n',2d0,1d0,
$ pars%g_eriov(isp),1,pars%npoles(isp),1,nri,
$ pars%g_eriov(jsp),1,nri,1,pars%npoles(jsp),
$ pars%g_apb,1+off1,pars%npoles(isp)+off1,
$ 1+off2,pars%npoles(jsp)+off2)
enddo
enddo
endif
call ga_add_diagonal(pars%g_apb,g_diag)
call ga_add_diagonal(g_amb,g_diag)
! Destroy ERIs
if (.not.ga_destroy(pars%g_eriov(1)))
& call errquit(pname//'could not destroy eriov GA',275,GA_ERR)
if (.not.ga_destroy(pars%g_eriov(2)))
& call errquit(pname//'could not destroy eriov GA',277,GA_ERR)
! First step is to get (A-B)^(1/2)
call ga_sync()
#ifdef SCALAPACK
info= ga_cholesky('L',g_amb)
#elif defined(PEIGS)
call ga_chol(g_amb, g_amb, info)
#else
call ga_chol_seq(g_amb, g_amb, info)
#endif
if (info.ne.0) call errquit("Cholesky Failed",298,GA_ERR)
! Transform (A+B) matrix with Cholesky factors
if(.not.ga_create(mt_dbl,npoles,npoles,'tmp',0,npoles,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('n','n',npoles,npoles,npoles,1d0,
& pars%g_apb,g_amb,0d0,g_tmp)
call ga_dgemm('t','n',npoles,npoles,npoles,1d0,
& g_amb,g_tmp,0d0,pars%g_apb)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
! Diagonalize (A+B) matrix
#ifdef PARALLEL_DIAG
#ifdef SCALAPACK
call dft_scaleig(pars%g_apb,pars%g_apb,omega,0)
#else
call ga_diag_std(pars%g_apb,pars%g_apb,omega)
#endif
#else
call ga_diag_std_seq(pars%g_apb,pars%g_apb,omega)
#endif
!
! Since we transformed the eigenvalue equation to an Hermitian
! form, we obtain the square of the neutral excitations
! (Equation (32)).
!
!$omp parallel do
do ipole=1,npoles
omega(ipole) = dsqrt(omega(ipole))
enddo
!$omp end parallel do
! Transform eigenvectors with (A-B)^(1/2)
if(.not.ga_create(mt_dbl,npoles,npoles,'tmp',0,npoles,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('n','n',npoles,npoles,npoles,1d0,
& g_amb,pars%g_apb,0d0,g_tmp)
! Gather neutral excitations vectors
call nga_distribution(g_diag,pars%me,ilo,ihi)
if(ihi.ge.ilo) then
call nga_put(g_diag,ilo,ihi,1d0/dsqrt(omega(ilo:ihi)),
& ihi-ilo+1)
endif
call ga_sync()
call ga_scale_cols(g_tmp,g_diag)
call ga_copy(g_tmp,pars%g_apb)
! Destroy unused GAs
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',344,GA_ERR)
if (.not.ga_destroy(g_amb))
& call errquit(pname//'could not destroy temp GA',346,GA_ERR)
if (.not.ga_destroy(g_diag))
& call errquit(pname//'could not destroy temp GA',348,GA_ERR)
endif
end subroutine

View file

@ -0,0 +1,200 @@
subroutine bse_analytic_tda(pars,wia,omega,maxpoles,totpoles,
$ singlet)
implicit none
#include "mafdecls.fh"
#include "errquit.fh"
#include "global.fh"
#include "bse.fh"
#ifdef SCALAPACK
#include "scaleig.fh"
integer ga_cholesky,ga_llt_i
external ga_cholesky,ga_llt_i
#endif
type(bse_params_t) :: pars
integer maxpoles, totpoles
double precision wia(maxpoles,pars%ipol)
double precision omega(totpoles)
logical singlet
character(*),parameter :: pname = 'bse_analytic_tda: '
integer nri,isp,jsp,npoles
integer ilo,ihi,jlo,jhi
integer nocc2(2),nvir2(2)
integer ipole,imo,jmo,amo,bmo
integer k,l,kproc,lproc,klocal,llocal,kglobal,lglobal
integer g_diag, g_tmp, g_amb, ktmp, ld, info, i, j
integer modtn,tovern,off1,off2
double precision alpha
integer ij,jb,llia,ulia
#ifdef USE_OPENMP
call util_blas_set_num_threads(pars%iMaxthreads)
#endif
npoles = pars%npoles(1)
if (pars%ipol.gt.1) npoles = npoles + pars%npoles(2)
nri = pars%nri
if(.not.ga_create(mt_dbl,totpoles,totpoles,'apb',0,npoles,
& pars%g_apb))
& call errquit(pname//'could not create A+B GA',0,GA_ERR)
if (.not.nga_create(mt_dbl,1,npoles,'wia',0,g_diag))
$ call errquit(pname//'could not create wia',0,GA_ERR)
if (pars%ipol.eq.1) then
nocc2(1) = pars%nocc(1)**2
nvir2(1) = pars%nvir(1)**2
call ga_distribution(pars%g_apb,pars%me,ilo,ihi,jlo,jhi)
if (jhi.ge.jlo) call nga_put(g_diag,ilo,ihi,wia,ihi-ilo+1)
! Screened Coulomb contribution
if(.not.ga_create(mt_dbl,nvir2(1),nocc2(1),'tmp',nvir2(1),0,
& g_tmp))
& call errquit(pname//'could not create temp GA',1,GA_ERR)
call ga_dgemm('t','n',nvir2(1),nocc2(1),pars%nri,-1d0,
% pars%g_erivv(1),pars%g_erioo(1),0d0,g_tmp)
call ga_distribution(g_tmp,pars%me,ilo,ihi,jlo,jhi)
call ga_access(g_tmp,ilo,ihi,jlo,jhi,ktmp,ld)
! Need to take care about indexing
do ij=jlo,jhi
jmo = (ij-1)/pars%nocc(1) + 1
imo = ij - (jmo-1)*pars%nocc(1)
do bmo = 1,pars%nvir(1)
llia = (imo-1)*pars%nvir(1) + 1
ulia = imo*pars%nvir(1)
jb = (jmo-1)*pars%nvir(1) + bmo
call ga_put(pars%g_apb,llia,ulia,jb,jb,
& dbl_mb(ktmp+(ij-jlo)*nvir2(1)+(bmo-1)*pars%nvir(1)),
& pars%nvir(1))
enddo
enddo
call ga_release(g_tmp,ilo,ihi,jlo,jhi)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
! Hartree contribution
if (singlet) then
call ga_dgemm('t','n',npoles,npoles,pars%nri,2d0,
& pars%g_eriov(1),pars%g_eriov(1),1d0,pars%g_apb)
endif
! Diagonal term
call ga_add_diagonal(pars%g_apb,g_diag)
! Destroy ERIs
if (.not.ga_destroy(pars%g_erioo(1)))
& call errquit(pname//'could not destroy erioo GA',1,GA_ERR)
if (.not.ga_destroy(pars%g_erivv(1)))
& call errquit(pname//'could not destroy erivv GA',1,GA_ERR)
if (.not.ga_destroy(pars%g_eriov(1)))
& call errquit(pname//'could not destroy eriov GA',1,GA_ERR)
! Diagonalize (A+B) matrix
#ifdef PARALLEL_DIAG
#ifdef SCALAPACK
call dft_scaleig(pars%g_apb,pars%g_apb,omega,0)
#else
call ga_diag_std(pars%g_apb,pars%g_apb,omega)
#endif
#else
call ga_diag_std_seq(pars%g_apb,pars%g_apb,omega)
#endif
! Destroy unused GAs
if (.not.ga_destroy(g_diag))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
else
nocc2(1) = pars%nocc(1)**2
nvir2(1) = pars%nvir(1)**2
nocc2(2) = pars%nocc(2)**2
nvir2(2) = pars%nvir(2)**2
! Gather eigenvalue differences
call ga_distribution(pars%g_eriov(1),pars%me,ilo,ihi,jlo,jhi)
call nga_put(g_diag,jlo,jhi,wia(1,1),jhi-jlo+1)
call ga_distribution(pars%g_eriov(2),pars%me,ilo,ihi,jlo,jhi)
call nga_put(g_diag,jlo+pars%npoles(1),jhi+pars%npoles(1),
$ wia(1,2),jhi-jlo+1)
! Form RPA matrix
if (.not.singlet) then
call ga_zero(pars%g_apb)
else
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
do jsp=1,pars%ipol
off2 = (jsp-1)*pars%npoles(1)
call ga_matmul_patch('t','n',1d0,0d0,
$ pars%g_eriov(isp),1,pars%npoles(isp),1,nri,
$ pars%g_eriov(jsp),1,nri,1,pars%npoles(jsp),
$ pars%g_apb,1+off1,pars%npoles(isp)+off1,
$ 1+off2,pars%npoles(jsp)+off2)
enddo
enddo
endif
call ga_add_diagonal(pars%g_apb,g_diag)
!Add Screened Coulomb contribution
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
if (.not.ga_destroy(pars%g_eriov(isp)))
& call errquit(pname//'could not destroy eriov GA',147,GA_ERR)
if (.not.ga_create(mt_dbl,nvir2(isp),nocc2(isp),'tmp',
$ nvir2(isp),0,g_tmp))
$ call errquit('bse_analytic_tda: could not create tmp',145,
$ GA_ERR)
call ga_dgemm('t','n',nvir2(isp),nocc2(isp),nri,-1d0,
$ pars%g_erivv(isp),pars%g_erioo(isp),0d0,g_tmp)
call ga_distribution(g_tmp,pars%me,ilo,ihi,jlo,jhi)
call ga_access(g_tmp,ilo,ihi,jlo,jhi,ktmp,ihi-ilo+1)
do ij=jlo,jhi
jmo = (ij-1)/pars%nocc(isp) + 1
imo = ij - (jmo-1)*pars%nocc(isp)
do bmo = 1,pars%nvir(isp)
llia = (imo-1)*pars%nvir(isp) + 1 + off1
ulia = imo*pars%nvir(isp) + off1
jb = (jmo-1)*pars%nvir(isp) + bmo + off1
call ga_acc(pars%g_apb,llia,ulia,jb,jb,
& dbl_mb(ktmp+(ij-jlo)*nvir2(isp)+
& (bmo-1)*pars%nvir(isp)),pars%nvir(isp),1d0)
enddo
enddo
call ga_release(g_tmp,ilo,ihi,jlo,jhi)
if (.not.ga_destroy(g_tmp))
& call errquit(pname//'could not destroy temp GA',172,GA_ERR)
! Destroy ERIs
if (.not.ga_destroy(pars%g_erioo(isp)))
& call errquit(pname//'could not destroy erioo GA',175,GA_ERR)
if (.not.ga_destroy(pars%g_erivv(isp)))
& call errquit(pname//'could not destroy erivv GA',177,GA_ERR)
enddo
! Diagonalize (A+B) matrix
#ifdef PARALLEL_DIAG
#ifdef SCALAPACK
call dft_scaleig(pars%g_apb,pars%g_apb,omega,0)
#else
call ga_diag_std(pars%g_apb,pars%g_apb,omega)
#endif
#else
call ga_diag_std_seq(pars%g_apb,pars%g_apb,omega)
#endif
! Destroy unused GAs
if (.not.ga_destroy(g_diag))
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
endif
end subroutine

78
src/bsemol/bse_buildw.F Normal file
View file

@ -0,0 +1,78 @@
subroutine bse_buildW(pars,wia,w,nmo,nri,maxpoles,ipol)
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "stdio.fh"
#include "bse.fh"
type(bse_params_t) :: pars
integer nmo, nri, maxpoles, ipol
character(*), parameter :: pname = 'bse_buildw: '
double precision wia(maxpoles,ipol)
double precision w(nri,nri)
integer :: ovlo,ovhi,kov(2),isp
integer :: ldpi, sizepi, iri, jri
integer :: alloc, info, ld, ipole
double precision,allocatable :: factor(:,:)
allocate(factor(maxpoles,ipol),stat=alloc)
if (alloc.ne.0)
& call errquit(pname//'allocation failed',0,MA_ERR)
C
C Decide the leading dimension of the \Pi array according to
C the rectangular full-packed (RFP) format
C
if (mod(nri,2).eq.0) then
ldpi = nri + 1
else
ldpi = nri
endif
sizepi = (nri*(nri+1))/2
! Get pointer to local integrals
do isp=1,pars%ipol
kov(isp) = 1
if (pars%mynpoles(isp).lt.1) cycle
ovlo = pars%ovlo(isp)
ovhi = pars%ovhi(isp)
call ga_access(pars%g_eriov(isp),1,nri,ovlo,ovhi,kov(isp),ld)
!$omp parallel do
do ipole=1,pars%mynpoles(isp)
factor(ipole,isp) = dsqrt(1.0d0/wia(ipole,isp))
enddo
!$omp end parallel do
enddo
! Build the Polarizability matrix and transform it to the
! dielectric matrix \epsilon = 1 - \pi
call gw_cdgw_buildpi('w',dbl_mb(kov(1)),dbl_mb(kov(2)),w,factor,
$ ldpi,nri,pars%mynpoles,pars%me.eq.0,
$ ipol,maxpoles)
call ga_dgop(1038,w,sizepi,'+')
!Factorize dieletric matrix for further use
call dpftrf('n','l',nri,w,info)
if(info.ne.0)
& call errquit(pname//'Failed to decompose dielectric',info,0)
!Deallocate
deallocate(factor,stat=alloc)
if (alloc.ne.0)
& call errquit(pname//'deallocation failed',0,MA_ERR)
!Release GAs
do isp=1,ipol
if (pars%mynpoles(isp).lt.1) cycle
ovlo = pars%ovlo(isp)
ovhi = pars%ovhi(isp)
call ga_release(pars%g_eriov(isp),1,nri,ovlo,ovhi)
enddo
end

569
src/bsemol/bse_davidson.F Normal file
View file

@ -0,0 +1,569 @@
subroutine bse_davidson(pars)
implicit none
#include "global.fh"
#include "errquit.fh"
#include "mafdecls.fh"
#include "bse.fh"
#include "cdft.fh"
type(bse_params_t) :: pars
character(*), parameter :: pname = 'bse_analytic: '
integer lSize,isp
integer l_omega, k_omega
integer l_dia,k_dia
integer lW,kW
integer maxpoles,maxnpoles,totnpoles,nri
integer ntrials,g_trials,npoles,ndim
integer g_x,g_y,g_ym,g_s,g_ks,g_mks,g_w,g_kw,g_mkw,off1
integer g_wia,g_tmp,info,lwork,kwork,idum,mdim,k_trials
integer k_x,k_y,k_ym,k_s,k_ks,k_mks,k_w,k_kw,k_mkw,k_r,l_r
integer l_map,k_map,l_tmp1,k_tmp1,l_tmp2,k_tmp2,k_x2,k_trials2
integer k_w2,l_tmp3,k_tmp3,l_tmp4,k_tmp4,k_kw2,k_mkw2
integer myld,vkv,l_vkv,myhi,mylo,l_ksw,ksw,jdum,myld2
integer subscript(2)
double precision minimum,norm(pars%nroots),rdum,factor
integer inode,nprocs,nrestart
integer ipole,jpole,itrial,imo,amo,idavidson
logical found,dohartree,stat,debug
double precision,external :: dasum
!debug = pars%me.eq.0
debug = .false.
! synchronize parallel processes
call ga_sync()
!initialize useful variables
nri = pars%nri
maxpoles = maxval(pars%mynpoles(1:pars%ipol))
maxnpoles = maxval(pars%npoles(1:pars%ipol))
totnpoles = maxnpoles*pars%ipol
npoles = sum(pars%npoles(1:pars%ipol))
ntrials = pars%nroots*pars%nspace
nprocs = ga_nnodes()
! allocation
if(.not.ma_push_get(mt_dbl,totnpoles,'dia',l_dia,k_dia))
& call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,npoles,'omega',l_omega,k_omega))
& call errquit(pname//'failed to allocate omega',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,npoles,'residual',l_r,k_r))
& call errquit(pname//'failed to allocate residual',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,maxpoles*pars%nmax,'tmp1',
& l_tmp1,k_tmp1))
& call errquit(pname//'failed to allocate tmp1',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,maxpoles*pars%nmax,'tmp2',
& l_tmp2,k_tmp2))
& call errquit(pname//'failed to allocate tmp2',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,maxpoles*pars%nmax,'tmp3',
& l_tmp3,k_tmp3))
& call errquit(pname//'failed to allocate tmp3',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,maxpoles*pars%nmax,'tmp4',
& l_tmp4,k_tmp4))
& call errquit(pname//'failed to allocate tmp4',0,MA_ERR)
if(.not.ma_push_get(mt_int,2*ga_nnodes()+1,'map',l_map,k_map))
& call errquit(pname//'failed to allocate map',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,ntrials*ntrials,'vkv',l_vkv,vkv))
& call errquit(pname//'failed to allocate vkv',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,nri**2,'W',lW,kW))
& call errquit(pname//'failed to allocate W',0,MA_ERR)
! Get map for irreg array
pars%nblocks = ga_nnodes()
do inode=0,ga_nnodes()-1
idum = pars%ovlo(1)
call ga_brdcst(mt_int,idum,ma_sizeof(mt_int,1,mt_byte),inode)
int_mb(k_map+inode) = idum
enddo
if (pars%ipol.gt.1) then
pars%nblocks = pars%nblocks + ga_nnodes()
do inode=0,ga_nnodes()-1
idum = pars%ovlo(2)
call ga_brdcst(mt_int,idum,ma_sizeof(mt_int,1,mt_byte),inode)
int_mb(k_map+ga_nnodes()+inode) = idum + pars%npoles(1)
enddo
endif
int_mb(k_map+pars%nblocks) = 1
! Allocate GA arrays
g_x = ga_create_handle()
call ga_set_data(g_x,2,(/npoles,ntrials/),mt_dbl)
call ga_set_array_name(g_x,'x')
call ga_set_tiled_irreg_proc_grid(g_x,int_mb(k_map),
$ (/pars%nblocks,1/),(/ga_nnodes(),1/))
stat = ga_allocate(g_x)
stat = stat.and.ga_duplicate(g_x,g_y,'y').and.
$ ga_duplicate(g_x,g_ym,'ym').and.
$ ga_duplicate(g_x,g_s,'s').and.
$ ga_duplicate(g_x,g_ks,'ks').and.
$ ga_duplicate(g_x,g_mks,'mks').and.
$ ga_duplicate(g_x,g_trials,'trials')
if(.not.stat) call errquit('ga create failed',103,GA_ERR)
call ga_zero(g_y)
call ga_zero(g_ym)
call ga_zero(g_s)
call ga_zero(g_ks)
call ga_zero(g_mks)
call ga_zero(g_trials)
! W, K*W, and M*K*W are fixed to nroots columns
g_w = ga_create_handle()
call ga_set_data(g_w, 2, (/npoles, pars%nroots/), mt_dbl)
call ga_set_array_name(g_w, 'w')
call ga_set_tiled_irreg_proc_grid(g_w, int_mb(k_map),
$ (/pars%nblocks, 1/), (/ga_nnodes(), 1/))
stat = ga_allocate(g_w)
stat = stat .and. ga_duplicate(g_w, g_kw, 'kw') .and.
$ ga_duplicate(g_w, g_mkw, 'mkw')
if(.not.stat) call errquit('Could not allocate g_w', 0, GA_ERR)
call ga_zero(g_w)
call ga_zero(g_kw)
call ga_zero(g_mkw)
g_wia = ga_create_handle()
call ga_set_data(g_wia,1,npoles,mt_dbl)
call ga_set_array_name(g_wia,'wia')
call ga_set_tiled_irreg_proc_grid(g_wia,int_mb(k_map),
$ pars%nblocks,ga_nnodes())
stat = stat.and.ga_allocate(g_wia)
if(.not.stat) call errquit('ga create failed',110,GA_ERR)
call ga_zero(g_wia)
! obtain "old" eigenvalue difference to build W
do isp=1,pars%ipol
lSize = (isp-1)*maxpoles
call bse_get_eia(dbl_mb(pars%k_mf_evals+(isp-1)*nbf_ao),
& dbl_mb(k_dia+lSize),pars%nocc(isp),
& pars%nvir(isp),pars%ovlo(isp),pars%ovhi(isp))
enddo
! obtain screened Coulomb matrix
call bse_buildw(pars,dbl_mb(k_dia),dbl_mb(kW),pars%nmo,nri,
& maxpoles,pars%ipol)
! transform ERIs using screened Coulomb matrix
call bse_wmn(pars,dbl_mb(kW),pars%nmo,nri,pars%ipol)
if(.not.ma_chop_stack(lW))
& call errquit(pname//'failed to chop stack',0,MA_ERR)
! obtain "GW" eigenvalue difference to build Casida matrix
do isp=1,pars%ipol
lSize = (isp-1)*maxpoles
call bse_get_eia(dbl_mb(pars%k_gw_evals+(isp-1)*nbf_ao),
& dbl_mb(k_dia+lSize),pars%nocc(isp),
& pars%nvir(isp),pars%ovlo(isp),
& pars%ovhi(isp))
enddo
! Get diagonal of Casida matrix
call bse_davidson_diagonal(pars,dbl_mb(k_dia),maxpoles)
! Get Eigenvalue differences
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
lSize = (isp-1)*maxpoles
call nga_put(g_wia,pars%ovlo(isp)+off1,pars%ovhi(isp)+off1,
$ dbl_mb(k_dia+lSize),pars%ovhi(isp)-pars%ovlo(isp)+1)
enddo
! Get guess vectos
if(pars%me.eq.0)
$ call bse_davidson_guess(pars,dbl_mb(k_dia),g_trials,npoles,
$ ntrials,pars%nroots)
call ga_brdcst(mt_dbl,dbl_mb(k_dia),
$ ma_sizeof(mt_dbl,npoles,mt_byte),0)
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
call dcopy(pars%mynpoles(isp),
$ dbl_mb(k_dia+off1+pars%ovlo(isp)-1),
$ 1,dbl_mb(k_dia+(isp-1)*pars%mynpoles(1)),1)
enddo
do isp=1,pars%ipol
if(.not.ga_create(mt_dbl,nri*pars%nvir(isp),pars%nvir(isp),
$ 'g_erivv',0,pars%nvir(isp),g_tmp))
$ call errquit('could not create GA',0,GA_ERR)
call ga_zero(g_tmp)
call ga_copy_patch('n',pars%g_erivv(isp),1,nri,1,
$ pars%nvir(isp)*pars%nvir(isp),
$ g_tmp,1,nri*pars%nvir(isp),1,
$ pars%nvir(isp))
if(.not.ga_destroy(pars%g_erivv(isp)))
$ call errquit('could not destroy GA',0,GA_ERR)
pars%g_erivv(isp) = g_tmp
enddo
! Start K-Davidson Algorithm
idavidson = 0
! Restarts should land here
1 continue
! K-orthogonalization
call bse_davidson_kvec(pars,g_trials,g_wia,g_x,npoles,
$ ntrials,pars%singlet,1d0)
call nga_access_block_grid(g_x, (/pars%me,0/), k_x, myld)
call nga_access_block_grid(g_trials, (/pars%me,0/),k_trials,myld)
if (pars%ipol.eq.1) then
myld2 = 0
call dcopy(myld*ntrials,dbl_mb(k_x),1,dbl_mb(k_tmp1),1)
call dcopy(myld*ntrials,dbl_mb(k_trials),1,dbl_mb(k_tmp2),1)
else
call nga_access_block_grid(g_x,(/pars%me+nprocs,0/),k_x2,myld2)
call nga_access_block_grid(g_trials,(/pars%me+nprocs,0/),
$ k_trials2,myld2)
do itrial=1,ntrials
call dcopy(myld,dbl_mb(k_x+(itrial-1)*myld),1,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_x2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)+myld),1)
call dcopy(myld,dbl_mb(k_trials+(itrial-1)*myld),1,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_trials2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)+myld),1)
enddo
endif
call bse_davidson_kortho(pars,dbl_mb(k_tmp2),dbl_mb(k_tmp1),
$ dbl_mb(k_tmp1),dbl_mb(vkv),dbl_mb(k_tmp1),
$ myld+myld2,ntrials,mdim,.true.)
if (pars%ipol.eq.1) then
call dcopy(myld*ntrials,dbl_mb(k_tmp2),1,dbl_mb(k_trials),1)
else
do itrial=1,ntrials
call dcopy(myld,dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)),1,
$ dbl_mb(k_trials+(itrial-1)*myld),1)
call dcopy(myld2,dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)+myld),
$ 1,dbl_mb(k_trials2+(itrial-1)*myld2),1)
enddo
call nga_release_update_block_grid(g_trials,
$ (/pars%me+nprocs,0/))
call nga_release_block_grid(g_x,(/pars%me+nprocs,0/))
endif
call nga_release_update_block_grid(g_trials,(/pars%me,0/))
call nga_release_block_grid(g_x,(/pars%me,0/))
! Starting vectors
call ga_copy(g_trials,g_x)
call bse_davidson_kvec(pars,g_x,g_wia,g_y,npoles,
$ ntrials,pars%singlet,1d0)
call bse_davidson_kvec(pars,g_y,g_wia,g_ym,npoles,
$ ntrials,.false.,-1d0)
call ga_copy(g_x,g_s)
call ga_copy(g_y,g_ks)
call ga_copy(g_ym,g_mks)
if (pars%me.eq.0 .and. idavidson.eq.0) then
write(*,*)
write(*,*) "--------------------------------"
write(*,*) " Iter Ntrials Nconv Res "
write(*,*) "--------------------------------"
endif
2 idavidson = idavidson + 1
if(idavidson.gt.pars%maxiter) goto 3
! Build projected matrix
call dcopy(ntrials*ntrials,0d0,0,dbl_mb(vkv),1)
do isp=1,pars%ipol
subscript = (/pars%me+(isp-1)*ga_nnodes(),0/)
call nga_access_block_grid(g_ks,subscript,k_ks,myld)
call nga_access_block_grid(g_mks,subscript,k_mks,myld)
call dgemm('t','n',ntrials,ntrials,myld,1d0,dbl_mb(k_ks),
$ myld,dbl_mb(k_mks),myld,1d0,dbl_mb(vkv),ntrials)
call nga_release_block_grid(g_ks,subscript)
call nga_release_block_grid(g_mks,subscript)
enddo
call ga_dgop((/194/),dbl_mb(vkv),ntrials**2,'+')
if (debug) call ma_print(dbl_mb(vkv),ntrials,ntrials,'Proj')
! Diagonalize projected matrix
call dsyev('v','l',ntrials,dbl_mb(vkv),ntrials,dbl_mb(k_omega),
$ rdum,-1,info)
if(.not.ma_push_get(mt_dbl,int(rdum),'work',lwork,kwork))
& call errquit(pname//'failed to allocate work',199,MA_ERR)
call dsyev('v','l',ntrials,dbl_mb(vkv),ntrials,dbl_mb(k_omega),
$ dbl_mb(kwork),int(rdum),info)
if(.not.ma_chop_stack(lwork))
& call errquit(pname//'failed to allocate work',192,MA_ERR)
if (info.ne.0)
$ call errquit(pname//'Diagonalization failed',206,0)
if (debug) call ma_print(dbl_mb(vkv),ntrials,ntrials,'Eigen')
if (debug) call ma_print(dbl_mb(k_omega),ntrials,1,'Omega')
! Project back vectors
do isp=1,pars%ipol
subscript = (/pars%me+(isp-1)*nprocs,0/)
call nga_access_block_grid(g_s,subscript,k_s,myld)
call nga_access_block_grid(g_ks,subscript,k_ks,myld)
call nga_access_block_grid(g_mks,subscript,k_mks,myld)
call nga_access_block_grid(g_x,subscript,k_x,myld)
call nga_access_block_grid(g_y,subscript,k_y,myld)
call nga_access_block_grid(g_ym,subscript,k_ym,myld)
call dgemm('n','n',myld,ntrials,ntrials,1d0,dbl_mb(k_s),myld,
$ dbl_mb(vkv),ntrials,0d0,dbl_mb(k_x),myld)
call dgemm('n','n',myld,ntrials,ntrials,1d0,dbl_mb(k_ks),myld,
$ dbl_mb(vkv),ntrials,0d0,dbl_mb(k_y),myld)
call dgemm('n','n',myld,ntrials,ntrials,1d0,dbl_mb(k_mks),
$ myld,dbl_mb(vkv),ntrials,0d0,dbl_mb(k_ym),myld)
call nga_release_block_grid(g_s,subscript)
call nga_release_block_grid(g_ks,subscript)
call nga_release_block_grid(g_mks,subscript)
call nga_release_update_block_grid(g_x,subscript)
call nga_release_update_block_grid(g_y,subscript)
call nga_release_update_block_grid(g_ym,subscript)
enddo
! Get residual
subscript = (/pars%me,0/)
call nga_access_block_grid(g_x,subscript,k_x,myld)
call nga_access_block_grid(g_ym,subscript,k_trials,myld)
if (pars%ipol.eq.1) then
myld2 = 0
call dcopy(myld*ntrials,dbl_mb(k_x),1,dbl_mb(k_tmp1),1)
call dcopy(myld*ntrials,dbl_mb(k_trials),1,dbl_mb(k_tmp2),1)
else
subscript = (/pars%me+nprocs,0/)
call nga_access_block_grid(g_x,subscript,k_x2,myld2)
call nga_access_block_grid(g_ym,subscript,k_trials2,myld2)
do itrial=1,pars%nroots
call dcopy(myld,dbl_mb(k_x+(itrial-1)*myld),1,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_x2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)+myld),1)
call dcopy(myld,dbl_mb(k_trials+(itrial-1)*myld),1,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_trials2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)+myld),1)
enddo
endif
call bse_davidson_residual(pars%nroots,ntrials,dbl_mb(k_omega),
$ dbl_mb(k_tmp1),dbl_mb(k_tmp2),dbl_mb(k_r),
$ dbl_mb(k_dia),dbl_mb(k_tmp3),myld+myld2,ndim,norm)
call nga_release_block_grid(g_x,(/pars%me,0/))
call nga_release_block_grid(g_ym,(/pars%me,0/))
call nga_access_block_grid(g_w,(/pars%me,0/),k_w,myld)
if(pars%ipol.eq.1) then
call dcopy(pars%nroots*myld,dbl_mb(k_tmp3),1,dbl_mb(k_w),1)
else
call nga_release_block_grid(g_x,subscript)
call nga_release_block_grid(g_ym,subscript)
call nga_access_block_grid(g_w,subscript,k_w2,myld2)
do itrial=1,pars%nroots
call dcopy(myld,dbl_mb(k_tmp3+(itrial-1)*(myld+myld2)),1,
$ dbl_mb(k_w+(itrial-1)*myld),1)
call dcopy(myld2,dbl_mb(k_tmp3+
$ (itrial-1)*(myld+myld2)+myld),1,
$ dbl_mb(k_w2+(itrial-1)*myld2),1)
enddo
call nga_release_update_block_grid(g_w,subscript)
endif
call nga_release_update_block_grid(g_w,(/pars%me,0/))
! Print info
if (pars%me.eq.0) then
write(*,9010) idavidson,ntrials,pars%nroots-ndim,maxval(norm)
endif
9010 format(I5,3X,I5,3X,I5,G16.5)
! All residuals are below the threshold
if (ndim.eq.0) goto 3
! Restart calculation if dimensions exceeded
if (ndim+ntrials.gt.pars%nmax) then
nrestart = min(ntrials, pars%nroots*pars%nspace)
call bse_davidson_restart(ntrials,nrestart,npoles,
$ g_x,g_y,g_ym,g_s,g_ks,g_mks,g_trials,l_vkv,vkv,
$ int_mb(k_map),pars%nblocks)
goto 1
endif
! Project-out current vectors
if(.not.ma_push_get(mt_dbl,ndim*ntrials,'ksw',l_ksw,ksw))
$ call errquit(pname//'ksw allocation failed',264,MA_ERR)
call dcopy(ntrials*ndim,0d0,0,dbl_mb(ksw),1)
do isp=1,pars%ipol
subscript = (/pars%me+(isp-1)*nprocs,0/)
call nga_access_block_grid(g_ks,subscript,k_ks,myld)
call nga_access_block_grid(g_w,subscript,k_w,myld)
call dgemm('t','n',ntrials,ndim,myld,1d0,dbl_mb(k_ks),myld,
$ dbl_mb(k_w),myld,1d0,dbl_mb(ksw),ntrials)
call nga_release_block_grid(g_ks,subscript)
call nga_release_block_grid(g_w,subscript)
enddo
call ga_dgop((/268/),dbl_mb(ksw),ntrials*ndim,'+')
do isp=1,pars%ipol
subscript = (/pars%me+(isp-1)*nprocs,0/)
call nga_access_block_grid(g_s,subscript,k_s,myld)
call nga_access_block_grid(g_w,subscript,k_w,myld)
call dgemm('n','n',myld,ndim,ntrials,-1d0,dbl_mb(k_s),myld,
$ dbl_mb(ksw),ntrials,1d0,dbl_mb(k_w),myld)
call nga_release_update_block_grid(k_w,subscript)
call nga_release_block_grid(k_s,subscript)
enddo
if(.not.ma_chop_stack(l_ksw))
$ call errquit(pname//'ksw deallocation failed',272,MA_ERR)
! K-orthogonalization
call bse_davidson_kvec(pars,g_w,g_wia,g_kw,npoles,
$ ndim,pars%singlet,1d0)
call bse_davidson_kvec(pars,g_kw,g_wia,g_mkw,npoles,
$ ndim,.false.,-1d0)
subscript = (/pars%me,0/)
call nga_access_block_grid(g_w,subscript,k_w,myld)
call nga_access_block_grid(g_kw,subscript,k_kw,myld)
call nga_access_block_grid(g_mkw,subscript,k_mkw,myld)
if (pars%ipol.eq.1) then
myld2 = 0
call dcopy(ndim*myld,dbl_mb(k_w),1,dbl_mb(k_tmp1),1)
call dcopy(ndim*myld,dbl_mb(k_kw),1,dbl_mb(k_tmp2),1)
call dcopy(ndim*myld,dbl_mb(k_mkw),1,dbl_mb(k_tmp3),1)
else
subscript = (/pars%me+nprocs,0/)
call nga_access_block_grid(g_w,subscript,k_w2,myld2)
call nga_access_block_grid(g_kw,subscript,k_kw2,myld2)
call nga_access_block_grid(g_mkw,subscript,k_mkw2,myld2)
do itrial=1,ndim
call dcopy(myld,dbl_mb(k_w+(itrial-1)*myld),1,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_w2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)+myld),1)
call dcopy(myld,dbl_mb(k_kw+(itrial-1)*myld),1,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_kw2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)+myld),1)
call dcopy(myld,dbl_mb(k_mkw+(itrial-1)*myld),1,
$ dbl_mb(k_tmp3+(itrial-1)*(myld+myld2)),1)
call dcopy(myld2,dbl_mb(k_mkw2+(itrial-1)*myld2),1,
$ dbl_mb(k_tmp3+(itrial-1)*(myld+myld2)+myld),1)
enddo
endif
call bse_davidson_kortho(pars,dbl_mb(k_tmp1),dbl_mb(k_tmp2),
$ dbl_mb(k_tmp3),dbl_mb(vkv),dbl_mb(k_tmp4),myld+myld2,
$ ndim,mdim,.false.)
if(pars%ipol.eq.1) then
call dcopy(mdim*myld,dbl_mb(k_tmp1),1,dbl_mb(k_w),1)
call dcopy(mdim*myld,dbl_mb(k_tmp2),1,dbl_mb(k_kw),1)
call dcopy(mdim*myld,dbl_mb(k_tmp3),1,dbl_mb(k_mkw),1)
else
do itrial=1,mdim
call dcopy(myld,dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)),1,
$ dbl_mb(k_w+(itrial-1)*myld),1)
call dcopy(myld2,
$ dbl_mb(k_tmp1+(itrial-1)*(myld+myld2)+myld),1,
$ dbl_mb(k_w2+(itrial-1)*myld2),1)
call dcopy(myld,dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)),1,
$ dbl_mb(k_kw+(itrial-1)*myld),1)
call dcopy(myld2,
$ dbl_mb(k_tmp2+(itrial-1)*(myld+myld2)+myld),1,
$ dbl_mb(k_kw2+(itrial-1)*myld2),1)
call dcopy(myld,dbl_mb(k_tmp3+(itrial-1)*(myld+myld2)),1,
$ dbl_mb(k_mkw+(itrial-1)*myld),1)
call dcopy(myld2,
$ dbl_mb(k_tmp3+(itrial-1)*(myld+myld2)+myld),1,
$ dbl_mb(k_mkw2+(itrial-1)*myld2),1)
enddo
call nga_release_update_block_grid(g_w,subscript)
call nga_release_update_block_grid(g_kw,subscript)
call nga_release_update_block_grid(g_mkw,subscript)
endif
call nga_release_update_block_grid(g_w,(/pars%me,0/))
call nga_release_update_block_grid(g_kw,(/pars%me,0/))
call nga_release_update_block_grid(g_mkw,(/pars%me,0/))
! All trial vectors were linearly-dependent
! restart calculation with current best guess
if (mdim.lt.3) then
nrestart = min(ntrials, pars%nroots*pars%nspace)
call bse_davidson_restart(ntrials,nrestart,npoles,
$ g_x,g_y,g_ym,g_s,g_ks,g_mks,g_trials,l_vkv,vkv,
$ int_mb(k_map),pars%nblocks)
goto 1
endif
ndim = mdim
! Extend the size of the current vectors
call bse_davidson_newvecs(ntrials,ndim,npoles,g_x,g_y,g_ym,
$ g_s,g_ks,g_mks,g_w,g_kw,g_mkw,
$ l_vkv,vkv,int_mb(k_map),pars%nblocks)
goto 2
3 continue
if(pars%me.eq.0) write(*,*)
stat = ga_destroy(pars%g_apb) .and. ga_destroy(pars%g_amb)
if(.not.stat) call errquit('deallocation failed',481,GA_ERR)
stat = ga_create(mt_dbl,npoles,pars%nroots,'X',0,pars%nroots,
$ pars%g_apb)
stat = ga_create(mt_dbl,npoles,pars%nroots,'Y',0,pars%nroots,
$ pars%g_amb)
call ga_zero(pars%g_apb)
call ga_zero(pars%g_amb)
if(.not.stat) call errquit('failed to creare GA',484,GA_ERR)
call nga_copy_patch('n',g_x,(/1,1/),(/npoles,pars%nroots/),
$ pars%g_apb,(/1,1/),(/npoles,pars%nroots/))
call nga_copy_patch('n',g_y,(/1,1/),(/npoles,pars%nroots/),
$ pars%g_amb,(/1,1/),(/npoles,pars%nroots/))
stat = ga_destroy(g_x).and.ga_destroy(g_y).and.ga_destroy(g_ym)
stat = stat.and.ga_destroy(g_s).and.ga_destroy(g_ks)
stat = stat.and.ga_destroy(g_mks).and.ga_destroy(g_w)
stat = stat.and.ga_destroy(g_kw).and.ga_destroy(g_mkw)
stat = stat.and.ga_destroy(g_trials).and.ga_destroy(g_wia)
do isp=1,pars%ipol
stat = stat.and.ga_destroy(pars%g_erioo(isp))
stat = stat.and.ga_destroy(pars%g_eriov(isp))
stat = stat.and.ga_destroy(pars%g_erivv(isp))
stat = stat.and.ga_destroy(pars%g_wov(isp))
enddo
if(.not.stat) call errquit('deallocation failed',496,GA_ERR)
if (.not.pars%tda) then
do ipole=0,pars%nroots-1
dbl_mb(k_omega+ipole) = dsqrt(dbl_mb(k_omega+ipole))
call ga_scale_patch(pars%g_apb,1,npoles,ipole+1,ipole+1,
$ dsqrt(dbl_mb(k_omega+ipole)))
call ga_scale_patch(pars%g_amb,1,npoles,ipole+1,ipole+1,
$ 1.0/dsqrt(dbl_mb(k_omega+ipole)))
enddo
endif
call ga_add(0.5d0,pars%g_apb,0.5d0,pars%g_amb,pars%g_apb)
call ga_add(1d0,pars%g_apb,-1d0,pars%g_amb,pars%g_amb)
call bse_oscstr(pars,dbl_mb(k_omega),npoles)
stat = ma_chop_stack(l_dia)
if(.not.stat) call errquit('deallocation failed',501,MA_ERR)
! Destroy GAs
if (.not.ga_destroy(pars%g_apb))
& call errquit(pname//'could not destroy temp GA',519,GA_ERR)
if (.not.ga_destroy(pars%g_amb))
& call errquit(pname//'could not destroy temp GA',519,GA_ERR)
end subroutine bse_davidson

View file

@ -0,0 +1,98 @@
subroutine bse_davidson_diagonal(pars,wia,maxpoles)
implicit none
#include "mafdecls.fh"
#include "errquit.fh"
#include "global.fh"
#include "bse.fh"
#include "util.fh"
type(bse_params_t) :: pars
integer maxpoles, npoles
double precision wia(maxpoles,pars%ipol)
character(*),parameter :: pname = 'bse_casida: '
integer isp,jsp,nri,off1,off2
integer ilo,ihi,jlo,jhi
integer nocc2(2),nvir2(2)
integer ipole,imo,jmo,amo,bmo
integer k,l,kproc,lproc,klocal,llocal,kglobal,lglobal
integer g_diag, g_tmp, g_amb, ktmp, ld, info, i, j
integer ij,jb,llia,ulia,ltmp2,ltmp3,ktmp2,ktmp3
double precision temp(maxpoles),factor
double precision, external :: ddot
#ifdef USE_OPENMP
call util_blas_set_num_threads(pars%iMaxthreads)
#endif
nri = pars%nri
npoles = pars%npoles(1)
if (pars%ipol.gt.1) npoles = npoles + pars%npoles(2)
if (.not.nga_create(mt_dbl,1,npoles,'apb',0,pars%g_apb))
$ call errquit(pname//'could not create wia',0,GA_ERR)
if (.not.nga_create(mt_dbl,1,npoles,'amb',0,pars%g_amb))
$ call errquit(pname//'could not create wia',0,GA_ERR)
if(.not.ma_push_get(mt_dbl,nri,'dia',ltmp2,ktmp2))
& call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR)
if(.not.ma_push_get(mt_dbl,nri,'dia',ltmp3,ktmp3))
& call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR)
factor = 2d0
if ((pars%singlet) .and. (pars%ipol.eq.1)) factor = 4d0
if ((pars%triplet) .and. (pars%ipol.eq.1)) factor = 0d0
if (pars%tda) factor = factor/2d0
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
! Eigenvalue differences
call ga_distribution(pars%g_eriov(isp),pars%me,ilo,ihi,jlo,jhi)
call dcopy(jhi-jlo+1,wia(1,isp),1,temp,1)
! Screened Coulomb from A
do ipole=jlo,jhi
imo = (ipole-1)/pars%nvir(isp) + 1
amo = ipole - (imo-1)*pars%nvir(isp)
imo = imo + (imo-1)*pars%nocc(isp)
amo = amo + (amo-1)*pars%nvir(isp)
call ga_get(pars%g_erioo(isp),1,nri,imo,imo,dbl_mb(ktmp2),nri)
call ga_get(pars%g_erivv(isp),1,nri,amo,amo,dbl_mb(ktmp3),nri)
temp(ipole-jlo+1) = temp(ipole-jlo+1) -
$ ddot(nri,dbl_mb(ktmp2),1,dbl_mb(ktmp3),1)
enddo
call nga_put(pars%g_apb,jlo+off1,jhi+off1,temp,jhi-jlo+1)
call nga_put(pars%g_amb,jlo+off1,jhi+off1,temp,jhi-jlo+1)
! Screened Coulomb from B
if (.not.pars%tda) then
call ga_access(pars%g_wov(isp),ilo,ihi,jlo,jhi,ktmp,ld)
do ipole=jlo,jhi
temp(ipole-jlo+1) = ddot(nri,dbl_mb(ktmp+(ipole-jlo)*nri),1,
$ dbl_mb(ktmp+(ipole-jlo)*nri),1)
enddo
call nga_acc(pars%g_apb,jlo+off1,jhi+off1,temp,jhi-jlo+1,-1d0)
call nga_acc(pars%g_amb,jlo+off1,jhi+off1,temp,jhi-jlo+1,1d0)
call ga_release(pars%g_wov(isp),ilo,ihi,jlo,jhi)
endif
! Hartree term
call ga_access(pars%g_eriov(isp),ilo,ihi,jlo,jhi,ktmp,ld)
do ipole=jlo,jhi
temp(ipole-jlo+1) =
$ factor*ddot(nri,dbl_mb(ktmp+(ipole-jlo)*nri),1,
$ dbl_mb(ktmp+(ipole-jlo)*nri),1)
enddo
call nga_acc(pars%g_apb,jlo+off1,jhi+off1,temp,jhi-jlo+1,1d0)
call ga_release(pars%g_eriov(isp),ilo,ihi,jlo,jhi)
enddo
if(.not.ma_chop_stack(ltmp2))
$ call errquit('bse_davidson_diagonal: chop stack',0,MA_ERR)
call ga_sync()
end subroutine

View file

@ -0,0 +1,141 @@
subroutine bse_davidson_guess(pars,diagonal,g_trials,npoles,
$ ntrials,nroots)
implicit none
#include "bse.fh"
#include "stdio.fh"
#include "mafdecls.fh"
#include "cdft.fh"
type(bse_params_t) :: pars
integer npoles,ntrials,nroots,g_trials
double precision diagonal(npoles)
double precision amb(npoles)
character*4 oname,vname
logical found
integer amo,imo,itrial,ipole,jpole,indeces(ntrials)
integer isym,jsym,i,j,pol,idum
integer allindeces(npoles)
double precision minimum
logical swapped
call nga_get(pars%g_apb,1,npoles,diagonal,npoles)
call nga_get(pars%g_amb,1,npoles,amb,npoles)
if (.not.pars%tda) then
diagonal(:) = diagonal(:)*amb(:)
endif
indeces(:) = -1
allindeces(:) = -1
do ipole=1,npoles
allindeces(ipole) = ipole
enddo
do ipole=npoles-1,1,-1
swapped = .false.
do jpole=1,ipole
if (diagonal(allindeces(jpole)).gt.
$ diagonal(allindeces(jpole+1))) then
idum = allindeces(jpole)
allindeces(jpole) = allindeces(jpole+1)
allindeces(jpole+1) = idum
swapped = .true.
endif
enddo
if(.not.swapped) exit
enddo
if (pars%dowindow) then
do ipole=1,npoles
if (dsqrt(diagonal(allindeces(ipole))).ge.pars%elower) then
if (ipole+ntrials-1.gt.npoles) then
indeces(1:ntrials) = allindeces(npoles-ntrials+1:npoles)
else
indeces(1:ntrials) = allindeces(ipole:ipole+ntrials-1)
endif
exit
endif
enddo
else
indeces(1:ntrials) = allindeces(1:ntrials)
endif
! Get the indeces of the smallest diagonal elements
! within energy window
#if 0
do itrial=1,ntrials
minimum = 1d6
do ipole=1,npoles
if (abs(diagonal(ipole)).lt.minimum) then
if (pars%dowindow) then
else
found = .true.
endif
do jpole=1,itrial
if (ipole.eq.indeces(jpole)) then
found = .false.
exit
endif
enddo
if (found) then
indeces(itrial) = ipole
minimum = abs(diagonal(ipole))
endif
endif
enddo
enddo
#endif
! Get symmetry of ground state
isym = 1
if (pars%ipol.eq.2) then
do i=1,pars%ipol
do j=1,pars%nocc(i)
isym = ieor(isym-1,int_mb(pars%k_irs(i)+j-1)-1)+1
enddo
enddo
endif
call sym_irrepname(geom, isym, oname)
! Start the guess vector
write(LuOut,*)
write(LuOut,9000) nroots
write(LuOut,9020)
write(LuOut,9010)
write(LuOut,9020)
do itrial=1,ntrials
ipole = indeces(itrial)
call ga_put(g_trials,ipole,ipole,itrial,itrial,1d0,1)
if(itrial.gt.nroots) cycle
! Print information
if (ipole.le.pars%npoles(1)) then
imo = (ipole-1)/pars%nvir(1) + 1
amo = ipole - (imo-1)*pars%nvir(1) + pars%nocc(1)
jsym = ieor(isym-1,int_mb(pars%k_irs(1)+imo-1)-1)+1
jsym = ieor(jsym-1,int_mb(pars%k_irs(1)+amo-1)-1)+1
pol = 1
else
imo = (ipole - pars%npoles(1) - 1)/pars%nvir(2) + 1
amo = ipole - pars%npoles(1) - (imo-1)*pars%nvir(2) +
$ pars%nocc(2)
jsym = ieor(isym-1,int_mb(pars%k_irs(2)+imo-1)-1)+1
jsym = ieor(jsym-1,int_mb(pars%k_irs(2)+amo-1)-1)+1
pol = 2
endif
call sym_irrepname(geom, jsym, vname)
write(luout,9030) itrial,pol,imo,amo,vname,
$ dsqrt(diagonal(ipole))*27.211396d0
enddo
9000 format(i5,' smallest diagonals (eV) ')
9010 format(' No. Spin Occ Vir Irrep Diag')
9020 format('------------------------------------')
9030 format(4i5,3x,a4,1x,f8.3)
end subroutine

View file

@ -0,0 +1,74 @@
subroutine bse_davidson_kortho(pars,v,kv,mkv,vkv,work,npoles,n,m,
$ first)
implicit none
#include "bse.fh"
#include "mafdecls.fh"
#include "errquit.fh"
type(bse_params_t) :: pars
logical first
integer npoles,n,m
integer l_work,k_work,info,i,j,k
double precision rdum
double precision v(npoles,n),kv(npoles,n),mkv(npoles,n)
double precision vkv(n,n),eig(n),work(npoles,n)
double precision sqrtvkv(n,n)
! V**T K V
call dgemm('t','n',n,n,npoles,1d0,v,npoles,kv,npoles,0d0,vkv,n)
call ga_dgop((/14/),vkv,n*n,'+')
! Allocate work field
call dsyev('v','l',n,vkv,n,eig,rdum,-1,info)
if(.not.ma_push_get(mt_dbl,int(rdum),'work',l_work,k_work))
$ call errquit('bse_kortho: failed to allocate work',199,MA_ERR)
! Diagonalize matrix
call dsyev('v','l',n,vkv,n,eig,dbl_mb(k_work),int(rdum),info)
if (info.ne.0)
$ call errquit("bse_davidson_kortho: diago failed",0,0)
! Deallocate field
if (.not.ma_chop_stack(l_work))
$ call errquit('bse_kortho: failed to deallocate work',199,MA_ERR)
! Remove small eigenpairs
m = n
k = 1
100 continue
do i=k,m
if (eig(i).lt.1d-9) then
do j=i+1,n
eig(j-1) = eig(j)
vkv(:,j-1) = vkv(:,j)
enddo
m = m - 1
k = i
goto 100
endif
enddo
! Obtain VKV^(-1/2)
do i=1,m
vkv(:,i) = vkv(:,i)/dsqrt(eig(i))
enddo
if (m.eq.0) return
! Obtain orthogonalized vectors
call dgemm('n','n',npoles,m,n,1d0,v,npoles,vkv,n,0d0,work,npoles)
call dcopy(npoles*m,work,1,v,1)
if (first) return
call dgemm('n','n',npoles,m,n,1d0,kv,npoles,vkv,n,0d0,work,npoles)
call dcopy(npoles*m,work,1,kv,1)
call dgemm('n','n',npoles,m,n,1d0,mkv,npoles,vkv,n,0d0,work,
$ npoles)
call dcopy(npoles*m,work,1,mkv,1)
end subroutine

View file

@ -0,0 +1,159 @@
subroutine bse_davidson_kvec(pars,vec,wia,kv,npoles,ntrials,
$ dohartree,sgn)
implicit none
#include "bse.fh"
#include "global.fh"
#include "mafdecls.fh"
#include "errquit.fh"
type(bse_params_t) :: pars
integer npoles,ntrials
integer vec,wia,kv,mv
double precision sgn
logical dohartree
integer ipol,npoles1,npoles2,uli,lli,ipole,off1
integer nocc(2),nvir(2),amo,imo,itrial,ula,lla,nri
integer ilo,ihi,jlo,jhi,ld,ld2, myld
integer l_tmp1,k_tmp1,k_eriov,isp,k_vec,k_erivv
integer g_tmp1, g_tmp2, g_tmp3, g_tmp4, k_tmp
integer subscript(2), k_tmp3
double precision temp(npoles),factor
ipol = pars%ipol
nri = pars%nri
nocc = pars%nocc
nvir = pars%nvir
npoles1 = nocc(1)*nvir(1)
npoles2 = nocc(2)*nvir(2)
call ga_zero(kv)
call ga_copy(vec,kv)
if (pars%tda .and. (sgn.eq.-1d0)) return
call ga_scale_rows(kv,wia)
! Screened Coulomb A
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
subscript = (/pars%me+(isp-1)*ga_nnodes(),0/)
if(.not.ga_create(mt_dbl,pars%nvir(isp),pars%nocc(isp),'g_tmp1',
$ 0,pars%nocc(isp),g_tmp1))
$ call errquit('bse_kvec: ga_create failed',42,GA_ERR)
call ga_zero(g_tmp1)
if(.not.ga_create(mt_dbl,nri,pars%nvir(isp),'g_tmp2',
$ nri,0,g_tmp2))
$ call errquit('bse_kvec: ga_create failed',45,GA_ERR)
call ga_zero(g_tmp2)
if(.not.ga_create(mt_dbl,nri*pars%nvir(isp),1,'g_tmp3',
$ 0,1,g_tmp3))
$ call errquit('bse_kvec: ga_create failed',48,GA_ERR)
call ga_zero(g_tmp3)
if(.not.ga_create(mt_dbl,pars%npoles(isp),1,'g_tmp4',
$ 0,1,g_tmp4))
$ call errquit('bse_kvec: ga_create failed',48,GA_ERR)
call ga_zero(g_tmp4)
call nga_access_block_grid(vec,subscript,k_vec,ld2)
call ga_distribution(pars%g_erivv(isp),pars%me,ilo,ihi,jlo,jhi)
call ga_access(pars%g_erivv(isp),ilo,ihi,jlo,jhi,k_erivv,ld)
call ga_access(g_tmp3,ilo,ihi,1,1,k_tmp3,ld)
do itrial=1,ntrials
call nga_copy_patch('n',vec,(/off1+1,itrial/),
$ (/off1+pars%npoles(isp),itrial/),g_tmp4,
$ (/1,1/),(/pars%npoles(isp),1/))
call ga_copy_patch('n',g_tmp4,1,pars%npoles(isp),1,1,
$ g_tmp1,1,pars%nvir(isp),1,pars%nocc(isp))
do imo=1,pars%nocc(isp)
lli = (imo-1)*pars%nocc(isp) + 1
uli = imo*pars%nocc(isp)
call ga_matmul_patch('n','t',1d0,0d0,
$ pars%g_erioo(isp),1,nri,lli,uli,
$ g_tmp1,1,pars%nocc(isp),1,pars%nvir(isp),
$ g_tmp2,1,nri,1,pars%nvir(isp))
call ga_copy_patch('n',g_tmp2,1,nri,1,pars%nvir(isp),
$ g_tmp3,1,nri*pars%nvir(isp),1,1)
call dgemv('t',ld,pars%nvir(isp),-1d0,dbl_mb(k_erivv),ld,
$ dbl_mb(k_tmp3),1,0d0,temp((imo-1)*pars%nvir(isp)+1),1)
enddo
call nga_acc(kv,(/1+off1,itrial/),
$ (/pars%npoles(isp)+off1,itrial/),temp,npoles,1d0)
enddo
if(.not. (ga_destroy(g_tmp1).and.ga_destroy(g_tmp2).and.
$ ga_destroy(g_tmp3).and.ga_destroy(g_tmp4)))
$ call errquit('bse_davidson_kvec: ga_destroy_failed',68,GA_ERR)
call ga_release(pars%g_erivv(isp),ilo,ihi,jlo,jhi)
call ga_release_update(g_tmp3,ilo,ihi,1,1)
call nga_release_block_grid(vec,subscript)
enddo
call ga_sync()
if (.not.ma_push_get(mt_dbl,nri*ntrials,'tmp1',l_tmp1,k_tmp1))
$ call errquit('bse_davidson_kvec: MA push failed',73,MA_ERR)
! Hartree
if (dohartree) then
factor = 4d0/pars%ipol
if (pars%tda) factor = factor/2d0
call dcopy(nri*ntrials,0d0,0,dbl_mb(k_tmp1),1)
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
subscript = (/pars%me+(isp-1)*ga_nnodes(),0/)
call ga_distribution(pars%g_eriov(isp),pars%me,ilo,ihi,jlo,jhi)
call ga_access(pars%g_eriov(isp),ilo,ihi,jlo,jhi,k_eriov,ld)
call nga_access_block_grid(vec,subscript,k_vec,ld)
call dgemm('n','n',nri,ntrials,ld,1d0,dbl_mb(k_eriov),nri,
$ dbl_mb(k_vec),ld,1d0,dbl_mb(k_tmp1),nri)
call ga_release(pars%g_eriov(isp),ilo,ihi,jlo,jhi)
call nga_release_block_grid(vec,subscript)
enddo
call ga_dgop((/93/),dbl_mb(k_tmp1),nri*ntrials,'+')
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
subscript = (/pars%me+(isp-1)*ga_nnodes(),0/)
call ga_distribution(pars%g_eriov(isp),pars%me,ilo,ihi,jlo,jhi)
call ga_access(pars%g_eriov(isp),ilo,ihi,jlo,jhi,k_eriov,ld)
call nga_access_block_grid(kv,subscript,k_vec,ld)
call dgemm('t','n',ld,ntrials,nri,factor,dbl_mb(k_eriov),nri,
$ dbl_mb(k_tmp1),nri,1d0,dbl_mb(k_vec),ld)
call ga_release(pars%g_eriov(isp),ilo,ihi,jlo,jhi)
call nga_release_update_block_grid(kv,subscript)
enddo
endif
! Screened Coulomb B
if(.not.pars%tda) then
do isp=1,pars%ipol
off1 = (isp-1)*pars%npoles(1)
subscript = (/pars%me+(isp-1)*ga_nnodes(),0/)
call ga_distribution(pars%g_wov(isp),pars%me,ilo,ihi,jlo,jhi)
call ga_access(pars%g_wov(isp),ilo,ihi,jlo,jhi,k_eriov,ld)
call nga_access_block_grid(vec,subscript,k_vec,ld)
call dgemm('n','n',nri,ntrials,ld,-1d0,dbl_mb(k_eriov),nri,
$ dbl_mb(k_vec),ld,0d0,dbl_mb(k_tmp1),nri)
call nga_release_block_grid(vec,subscript)
call ga_dgop((/115/),dbl_mb(k_tmp1),nri*ntrials,'+')
call nga_access_block_grid(kv,subscript,k_vec,ld)
call dgemm('t','n',ld,ntrials,nri,sgn,dbl_mb(k_eriov),nri,
$ dbl_mb(k_tmp1),nri,1d0,dbl_mb(k_vec),ld)
call ga_release(pars%g_wov(isp),ilo,ihi,jlo,jhi)
call nga_release_update_block_grid(kv,subscript)
enddo
endif
if (.not.ma_chop_stack(l_tmp1))
$ call errquit('bse_kvec: failed chop stack',122,MA_ERR)
end

View file

@ -0,0 +1,79 @@
subroutine bse_davidson_newvecs(ntrials,ndim,npoles,g_x,g_y,
$ g_ym,g_s,g_ks,g_mks,g_w,g_kw,g_mkw,lvkv,vkv,map,
$ nblocks)
implicit none
#include "global.fh"
#include "mafdecls.fh"
#include "errquit.fh"
integer ntrials,ndim,npoles,nblocks
integer g_x,g_y,g_ym,g_s,g_ks,g_mks
integer g_w,g_kw,g_mkw,lvkv,vkv
integer map(nblocks),grid(2)
logical stat
grid = (/ga_nnodes(),1/)
! Save working vectors
call ga_copy(g_s,g_x)
call ga_copy(g_ks,g_y)
call ga_copy(g_mks,g_ym)
! Destroy working vectors
stat = ga_destroy(g_s).and.ga_destroy(g_ks).and.ga_destroy(g_mks)
if(.not.stat)
$ call errquit('bse_davidson_newvecs: destroy failed',1,GA_ERR)
! Create new working vectors
g_s = ga_create_handle()
call ga_set_data(g_s,2,(/npoles,ntrials+ndim/),mt_dbl)
call ga_set_array_name(g_s,'s')
call ga_set_tiled_irreg_proc_grid(g_s,map,(/nblocks,1/),grid)
stat = ga_allocate(g_s)
stat = stat.and.ga_duplicate(g_s,g_ks,'ks').and.
$ ga_duplicate(g_s,g_mks,'mks')
if(.not.stat)
$ call errquit('bse_davidson_newvecs: create failed',1,GA_ERR)
! Copy back old working vectors plus new vectors
call ga_copy_patch('n',g_x,1,npoles,1,ntrials,
$ g_s,1,npoles,1,ntrials)
call ga_copy_patch('n',g_w,1,npoles,1,ndim,
$ g_s,1,npoles,ntrials+1,ntrials+ndim)
call ga_copy_patch('n',g_y,1,npoles,1,ntrials,
$ g_ks,1,npoles,1,ntrials)
call ga_copy_patch('n',g_kw,1,npoles,1,ndim,
$ g_ks,1,npoles,ntrials+1,ntrials+ndim)
call ga_copy_patch('n',g_ym,1,npoles,1,ntrials,
$ g_mks,1,npoles,1,ntrials)
call ga_copy_patch('n',g_mkw,1,npoles,1,ndim,
$ g_mks,1,npoles,ntrials+1,ntrials+ndim)
! Update ntrials
ntrials = ntrials + ndim
! Destroy old solutions
stat = ga_destroy(g_x).and.ga_destroy(g_y).and.ga_destroy(g_ym)
if(.not.stat)
$ call errquit('bse_davidson_newvecs: destroy failed',2,GA_ERR)
! Create new solutions
stat = ga_duplicate(g_s,g_x,'x').and.
$ ga_duplicate(g_s,g_y,'y').and.
$ ga_duplicate(g_s,g_ym,'ym')
if(.not.stat)
$ call errquit('bse_davidson_newvecs: create failed',2,GA_ERR)
! Reallocate VKV
if(.not.ma_chop_stack(lvkv))
$ call errquit('bse_davidson_restart: alloc failed',3,MA_ERR)
if(.not.ma_push_get(mt_dbl,ntrials*ntrials,'vkv',lvkv,vkv))
$ call errquit('bse_davidson_restart: dealloc failed',3,MA_ERR)
end subroutine

View file

@ -0,0 +1,32 @@
subroutine bse_davidson_residual(nroots,ntrials,omega,x,ym,r,dia,
$ w,npoles,ndim,norm)
implicit none
logical main
integer nroots,npoles,ntrials,ndim
double precision omega(nroots)
double precision ym(npoles,nroots),r(npoles)
double precision x(npoles,nroots),dia(npoles),w(npoles,nroots)
logical consecutive
integer iroot
double precision norm(nroots)
consecutive = .true.
ndim = 0
do iroot=1,nroots
r(:) = ym(:,iroot) - omega(iroot)*x(:,iroot)
norm(iroot) = maxval(abs(r))
call ga_dgop((/1/),norm(iroot),1,"maxabs")
if (consecutive .and. norm(iroot).gt.1d-4) consecutive = .false.
if ((norm(iroot).gt.1d-4) .or. (.not.consecutive)) then
ndim = ndim + 1
w(:,ndim) = r(:)/(dia(:) - omega(iroot))
endif
enddo
end subroutine

View file

@ -0,0 +1,78 @@
subroutine bse_davidson_restart(ntrials,nroots,npoles,g_x,g_y,
$ g_ym,g_s,g_ks,g_mks,g_trials,
$ lvkv,vkv,map,nblocks)
implicit none
#include "global.fh"
#include "mafdecls.fh"
#include "errquit.fh"
integer ntrials,nroots,npoles,nblocks
integer g_x,g_y,g_ym,g_s,g_ks,g_mks,lvkv,vkv
integer map(nblocks),grid(2),g_trials
logical stat
grid = (/ga_nnodes(),1/)
! Reset the number of trials
ntrials = nroots
! Copy current solutions
call ga_copy_patch('n',g_x,1,npoles,1,ntrials,
$ g_trials,1,npoles,1,ntrials)
! call ga_copy_patch('n',g_y,1,npoles,1,ntrials,
! $ g_ks,1,npoles,1,ntrials)
! call ga_copy_patch('n',g_ym,1,npoles,1,ntrials,
! $ g_mks,1,npoles,1,ntrials)
! Destroy current solutions
stat = ga_destroy(g_x).and.ga_destroy(g_y).and.ga_destroy(g_ym)
if (.not.stat)
$ call errquit('bse_davidson_restar: failed destroy',1,GA_ERR)
! Re-create new solutions
g_x = ga_create_handle()
call ga_set_data(g_x,2,(/npoles,ntrials/),mt_dbl)
call ga_set_array_name(g_x,'x')
call ga_set_tiled_irreg_proc_grid(g_x,map,(/nblocks,1/),grid)
stat = ga_allocate(g_x)
stat = stat.and.ga_duplicate(g_x,g_y,'y').and.
$ ga_duplicate(g_x,g_ym,'ym')
if (.not.stat)
$ call errquit('bse_davidson_restar: failed create',1,GA_ERR)
! Copy back solutions
! call ga_copy_patch('n',g_s,1,npoles,1,ntrials,
! $ g_x,1,npoles,1,ntrials)
! call ga_copy_patch('n',g_ks,1,npoles,1,ntrials,
! $ g_y,1,npoles,1,ntrials)
! call ga_copy_patch('n',g_mks,1,npoles,1,ntrials,
! $ g_ym,1,npoles,1,ntrials)
! Destroy working vectors
stat = ga_destroy(g_s).and.ga_destroy(g_ks).and.ga_destroy(g_mks)
if (.not.stat)
$ call errquit('bse_davidson_restar: failed destroy',2,GA_ERR)
! Re-create new solutions
stat = ga_duplicate(g_x,g_s,'s').and.
$ ga_duplicate(g_x,g_ks,'ks').and.
$ ga_duplicate(g_x,g_mks,'mks')
if (.not.stat)
$ call errquit('bse_davidson_restar: failed create',2,GA_ERR)
! Copy current solutions
! call ga_copy(g_x,g_s)
! call ga_copy(g_y,g_ks)
! call ga_copy(g_ym,g_mks)
! Reallocate VKV
if(.not.ma_chop_stack(lvkv))
$ call errquit('bse_davidson_restart: alloc failed',3,MA_ERR)
if(.not.ma_push_get(mt_dbl,ntrials*ntrials,'vkv',lvkv,vkv))
$ call errquit('bse_davidson_restart: dealloc failed',3,MA_ERR)
end subroutine

33
src/bsemol/bse_defaults.F Normal file
View file

@ -0,0 +1,33 @@
subroutine bse_defaults(rtdb)
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "rtdb.fh"
integer rtdb
logical ok
ok = rtdb_put(rtdb,'bse:analytic',mt_log,1,.true.) .and.
& rtdb_put(rtdb,'bse:davidson',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:lanczos',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:tda',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:ncap',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:ncapr',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:singlet',mt_log,1,.true.) .and.
& rtdb_put(rtdb,'bse:triplet',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:dowindow',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:elower',mt_dbl,1,0d0) .and.
& rtdb_put(rtdb,'bse:upper',mt_dbl,1,0d0) .and.
& rtdb_put(rtdb,'bse:nroots',mt_int,1,10) .and.
& rtdb_put(rtdb,'bse:nmax',mt_int,1,200) .and.
& rtdb_put(rtdb,'bse:skipgw',mt_log,1,.false.) .and.
& rtdb_put(rtdb,'bse:maxiter',mt_int,1,100) .and.
& rtdb_put(rtdb,'bse:nspace',mt_int,1,4) .and.
& rtdb_put(rtdb,'bse:init',mt_log,1,.true.)
if (.not.ok)
& call errquit('bse_defaults: could not put to rtdb',0,RTDB_ERR)
end subroutine bse_defaults

38
src/bsemol/bse_finalize.F Normal file
View file

@ -0,0 +1,38 @@
subroutine bse_finalize(pars)
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "bas.fh"
#include "cdft.fh"
#include "bse.fh"
type(bse_params_t) :: pars
character(*),parameter :: pname = 'bse_finalize: '
integer i
logical ok
do i=1,ipol
if (.not.ga_destroy(pars%g_movecs(i)))
$ call errquit(pname//'failed to destroy g_movecs',0,GA_ERR)
enddo
ok = ma_chop_stack(pars%l_mf_evals)
if (.not.ok)
& call errquit(pname//'could not chop stack',0,MA_ERR)
call schwarz_tidy()
call int_terminate()
call fock_2e_tidy(pars%rtdb)
if (.not.geom_destroy(geom))
$ call errquit(pname//'geom_destroy failed',0,GEOM_ERR)
if (.not.bas_destroy(ao_bas_han))
$ call errquit(pname//'bas_destroy failed',0,BASIS_ERR)
if (.not.bas_destroy(cd_bas_han))
$ call errquit(pname//'bas_destroy failed',0,BASIS_ERR)
end subroutine

15
src/bsemol/bse_get_eia.F Normal file
View file

@ -0,0 +1,15 @@
subroutine bse_get_eia(eig,eia,nocc,nvir,llpole,ulpole)
implicit none
integer,intent(in) :: nocc,nvir
double precision,intent(in),dimension(nocc+nvir) :: eig
double precision,intent(out),dimension(*) :: eia
integer ipole,llpole,ulpole,imo,amo
c
do ipole=llpole,ulpole
imo = (ipole-1)/nvir + 1
amo = ipole - (imo-1)*nvir + nocc
eia(ipole-llpole+1) = eig(amo) - eig(imo)
enddo
c
return
end subroutine

322
src/bsemol/bse_init.F Normal file
View file

@ -0,0 +1,322 @@
subroutine bse_init(pars)
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "cdft.fh"
#include "bas.fh"
#include "case.fh"
#include "geom.fh"
#include "global.fh"
#include "bse.fh"
#include "util.fh"
#include "stdio.fh"
#include "rtdb.fh"
#ifdef SCALAPACK
#include "scaleig.fh"
#endif
type(bse_params_t) :: pars
character(*), parameter :: pname = 'bse_init: '
character(len=255) :: basisname, scftype, gwevals
integer i, ilo, ihi, j, g_work
integer nbf_temp, ipol_temp, nocc, nvir, nmo_temp(2)
logical ldum, skipgw
logical int_normalize, int_norm_2c
external int_normalize, int_norm_2c
c
logical movecs_read, movecs_read_header, atom_tag_check
external movecs_read, movecs_read_header, atom_tag_check
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
integer resultlen,tid,cpu,ierr,iproc,itid
character(len=80) names
double precision t0, ehomo, vxddp, vxddm, aq2, zeta
integer,parameter :: unitno = 66
integer :: nmix(2)
logical ncapr
#ifdef USE_OPENMP
integer,external :: omp_get_max_threads, omp_get_thread_num
integer,external :: omp_get_num_threads
#endif
! Initialize GA variables
pars%me = ga_nodeid()
pars%nprocs = ga_nnodes()
! Initialize OMP variables
#ifdef USE_OPENMP
pars%iMaxThreads = omp_get_max_threads()
#else
pars%iMaxThreads = 1
#endif
pars%npoles(:) = 0
! Get defaults in case input block is not present
if (.not.rtdb_get(pars%rtdb,'bse:init',mt_log,1,ldum))
& call bse_defaults(pars%rtdb)
! Check method to solve BSE
if (.not.rtdb_get(pars%rtdb,'bse:analytic',mt_log,1,
& pars%analytic))
& call errquit(pname//'could not read from rtdb',1,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb,'bse:davidson',mt_log,1,
& pars%davidson))
& call errquit(pname//'could not read from rtdb',2,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb,'bse:lanczos',mt_log,1,
& pars%lanczos))
& call errquit(pname//'could not read from rtdb',3,RTDB_ERR)
if ((pars%analytic.and.(pars%davidson.or.pars%lanczos)) .or.
& (pars%davidson.and.(pars%analytic.or.pars%lanczos)) .or.
& (pars%lanczos.and.(pars%analytic.or.pars%davidson)))
& call errquit(pname//'more than one method',4,INPUT_ERR)
if (.not.rtdb_get(pars%rtdb, 'bse:skipgw', mt_log, 1, skipgw))
& call errquit(pname//'could not read skipgw',5, RTDB_ERR)
! Check wether we are computing singlet or triplet excitations
if(.not.rtdb_get(pars%rtdb,'bse:singlet',mt_log,1,pars%singlet))
& call errquit(pname//'could not read from rtdb',1,RTDB_ERR)
if(.not.rtdb_get(pars%rtdb,'bse:triplet',mt_log,1,pars%triplet))
& call errquit(pname//'could not read from rtdb',1,RTDB_ERR)
if((pars%singlet .and. pars%triplet).and.(pars%analytic)) then
call errquit(pname//'select either triplet or singlet',0,
& INPUT_ERR)
endif
! Check if we are doing Tamm-Dancoff approximation
if(.not.rtdb_get(pars%rtdb,'bse:tda',mt_log,1,pars%tda))
& call errquit(pname//'could not read from rtdb',101,RTDB_ERR)
! Check if Output is restricted to an energy window
if(.not.rtdb_get(pars%rtdb,'bse:dowindow',mt_log,1,pars%dowindow))
& call errquit(pname//'could not read from rtdb',105,RTDB_ERR)
if(pars%dowindow) then
if(.not.rtdb_get(pars%rtdb,'bse:elower',mt_dbl,1,pars%elower))
& call errquit(pname//'could not read from rtdb',108,RTDB_ERR)
if(.not.rtdb_get(pars%rtdb,'bse:eupper',mt_dbl,1,pars%eupper))
& call errquit(pname//'could not read from rtdb',110,RTDB_ERR)
endif
! Number of roots
if(.not.rtdb_get(pars%rtdb,'bse:nroots',mt_int,1,pars%nroots))
& call errquit(pname//'could not read from rtdb',115,RTDB_ERR)
! Starting space size
if(.not.rtdb_get(pars%rtdb, 'bse:nspace', mt_int, 1, pars%nspace))
& call errquit(pname//'could not read from rtdb', 119, RTDB_ERR)
! Maximum Number of roots
if(.not.rtdb_get(pars%rtdb,'bse:nmax',mt_int,1,pars%nmax))
& call errquit(pname//'could not read from rtdb',123,RTDB_ERR)
if ((pars%nmax.le.5*pars%nroots).and.pars%davidson) then
call errquit(pname//'maxroots should be at least 5*nroots',
& 126,INPUT_ERR)
endif
if((pars%nroots*pars%nspace.gt.pars%nmax).and.pars%davidson)
& call errquit(pname//'maxroots should be at least nroots*nspace',
& 130, INPUT_ERR)
! Maximum number of iterations
if(.not.rtdb_get(pars%rtdb,'bse:maxiter',mt_int,1,pars%maxiter))
& call errquit(pname//'could not read from rtdb',134,RTDB_ERR)
! Get some info from DFT module
call dft_rdinput(pars%rtdb)
if (cdfit) then
if (.not.bas_destroy(cd_bas_han))
$ call errquit(pname//'failed to destroy cd_bas',0,BASIS_ERR)
endif
if (.not.rtdb_get(pars%rtdb, 'dft:itol2e',mt_int,1,itol2e))
$ call errquit(pname//'failed to read itol2e',0,RTDB_ERR)
pars%tol2e = 10.0d0**(-itol2e-2)
! Fitting basis
if (.not.bas_create(cd_bas_han,'ri basis'))
$ call errquit(pname//'bas_create failed',0,BASIS_ERR)
if (.not.bas_rtdb_load(pars%rtdb,geom,cd_bas_han,'ri basis'))then
if (.not.bas_rtdb_load(pars%rtdb,geom,cd_bas_han,'cd basis'))
$ call errquit(pname//'a "cd basis" or "ri basis" is needed',0,
$ BASIS_ERR)
endif
call int_init(pars%rtdb, 1, cd_bas_han)
if (.not.int_norm_2c(pars%rtdb, cd_bas_han))
$ call errquit(pname//'int_norm_2c failed',0,INT_ERR)
call int_terminate()
if (.not. bas_numbf(cd_bas_han, nbf_cd))
$ call errquit(pname//'basis set error',0,BASIS_ERR)
if (.not. bas_nprim_cn_max(cd_bas_han,nbf_cd_mxprim))
$ call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_high_angular(cd_bas_han,nbf_cd_mxang))
$ call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_ncontr_cn_max(cd_bas_han,nbf_cd_mxcont))
$ call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_nbf_cn_max(cd_bas_han,nbf_cd_mxnbf_cn))
$ call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_nbf_ce_max(cd_bas_han,nbf_cd_mxnbf_ce))
$ call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_numcont(cd_bas_han,nshells_cd))
$ call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (bas_is_spherical(ao_bas_han).and.
$ (.not.bas_is_spherical(cd_bas_han)))
$ call int_app_set_no_texas(pars%rtdb)
pars%nri = nbf_cd
!Spin multiplicity, occupations
if (.not.rtdb_get(pars%rtdb,'dft:ipol',mt_int,1,ipol))
$ call errquit(pname//'failed to read ipol',0,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb, 'dft:noc',mt_int,2,noc))
$ call errquit(pname//'failed to read noc',0,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb,'dft:mult',mt_int,1,mult))
$ call errquit(pname//'failed to read mult',0,RTDB_ERR)
pars%ipol = ipol
! MO vectors
if (.not. rtdb_cget(pars%rtdb, 'dft:output vectors',1,movecs_out))
$ call errquit(pname//'no final MO found',0,RTDB_ERR)
if (.not.ma_push_get(mt_dbl,nbf_ao*ipol,
$ 'mf evals', pars%l_mf_evals, pars%k_mf_evals))
$ call errquit(pname//'failed to allocate evals',0,MA_ERR)
if (.not.ma_push_get(mt_dbl,nbf_ao*ipol,
$ 'gw evals', pars%l_gw_evals, pars%k_gw_evals))
$ call errquit(pname//'failed to allocate evals',0,MA_ERR)
if (.not.ma_push_get(mt_dbl,nbf_ao*ipol,'occupancies', pars%l_occ,
$ pars%k_occ))
$ call errquit(pname//'failed to allocate occ',0,MA_ERR)
if (.not.movecs_read_header(movecs_out,title,basisname,
$ scftype,nbf_temp,ipol_temp,nmo_temp,2))
$ call errquit(pname//'failed to read MO header',0,DISK_ERR)
if (nbf_ao.ne.nbf_temp)
$ call errquit(pname//'corrupted MO vectors',0,DISK_ERR)
if (ipol.ne.ipol_temp)
$ call errquit(pname//'corrupted MO vectors',0,DISK_ERR)
pars%nmo = nmo_temp(1)
do i=1,ipol
pars%g_movecs(i) = ga_create_atom_blocked(geom,
$ ao_bas_han,'MO eigenvectors')
if (.not.movecs_read(movecs_out,i,
$ dbl_mb(pars%k_occ+(i-1)*nbf_ao),
$ dbl_mb(pars%k_mf_evals+(i-1)*nbf_ao),
$ pars%g_movecs(i)))
$ call errquit(pname//'failed to read MO vectors',0,DISK_ERR)
nocc = 0
nvir = 0
do j=1,pars%nmo
if (dbl_mb(pars%k_occ+(i-1)*nbf_ao+j-1).gt.0.d0) then
nocc = nocc + 1
cycle
else
nvir = nvir + 1
endif
enddo
pars%nocc(i) = nocc
pars%nvir(i) = nvir
if (nocc+nvir.ne.pars%nmo)
$ call errquit(pname//'something went wrong',0,0)
pars%npoles(i) = pars%nocc(i)*pars%nvir(i)
enddo
if ((pars%nroots.eq.-1).or.(pars%nroots.gt.sum(pars%npoles))) then
pars%nroots = sum(pars%npoles)
endif
if (pars%nmax.gt.sum(pars%npoles)) then
pars%nmax = sum(pars%npoles)
endif
! Obtain NCAP DD shift
1001 format(2X,A12,F8.3,A3)
if (.not.rtdb_get(pars%rtdb,'bse:ncap',mt_log,1,pars%ncap))
& call errquit(pname//'rtdb_get failed',176,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb,'bse:ncapr',mt_log,1,ncapr))
& call errquit(pname//'rtdb_get failed',176,RTDB_ERR)
pars%ncap = ncapr .or. pars%ncap
if (pars%ncap) then
zeta = 0.304121d0
if (abs(xfac(87)).gt.1d-10) zeta=0.5d0
aq2 = 2d0/9d0 * (1d0-zeta)**2
do i=1,ipol
ehomo = dbl_mb(pars%k_mf_evals+(i-1)*nbf_ao+pars%nocc(i)-1)
vxddp = -0.5d0*aq2*(1d0 + sqrt(1d0-4d0*ehomo/aq2))!*0.85d0
vxddm = -0.5d0*aq2*(1d0 - sqrt(1d0-4d0*ehomo/aq2))!*0.85d0
pars%vxddp(i) = vxddp
pars%vxddm(i) = vxddm
if (pars%me.eq.0) then
write(luout,*) ' Shifting eigvals by 100% NCAP DD'
write(luout,1001) ' v_x^{DD-}: ',vxddp*ha2ev, ' eV'
write(luout,1001) ' v_x^(DD+): ',vxddm*ha2ev, ' eV'
endif
call daxpy(pars%nocc(i),1d0,vxddp,0,
& dbl_mb(pars%k_mf_evals+(i-1)*nbf_ao),1)
call daxpy(pars%nvir(i),1d0,vxddm,0,
& dbl_mb(pars%k_mf_evals+(i-1)*nbf_ao+pars%nocc(i)),1)
enddo
endif
! Read GW QP energies
!if (.not.pars%ncap) then
if(skipgw) then
call dcopy(nbf_ao*pars%ipol,dbl_mb(pars%k_mf_evals),1,
& dbl_mb(pars%k_gw_evals),1)
else
call util_file_name('gwevals',.false.,.false.,gwevals)
call util_file_name_resolve(gwevals,.false.)
if (pars%me.eq.0) then
open(unit=unitno,status='old',form='unformatted',file=gwevals)
do i=1,pars%ipol
call sread(unitno, dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao),
& pars%nmo)
enddo
close(unit=unitno)
endif
call ga_brdcst(mt_dbl,dbl_mb(pars%k_gw_evals),
& ma_sizeof(mt_dbl,nbf_ao*pars%ipol,mt_byte),0)
endif
!3-center ERIs
call int_init(pars%rtdb, 2, (/ao_bas_han, cd_bas_han/))
call print_integrals((/ao_bas_han,cd_bas_han/),.false.)
call schwarz_init(geom, ao_bas_han)
call scf_get_fock_param(pars%rtdb, pars%tol2e)
! Obtain irreps information
do i=1,ipol
if(.not.ma_push_get(mt_int,nbf_ao,'irreps',pars%l_irs(i),
$ pars%k_irs(i))) call errquit
$ ('bse_init: failed to allocate irs',0,MA_ERR)
if(.not.ga_create(mt_dbl,nbf_ao,pars%nmo,'Work',nbf_ao,-1,
$ g_work))
$ call errquit('bse_init: failed to allocate work',0,GA_ERR)
call ga_copy_patch('n',pars%g_movecs(i),1,nbf_ao,1,pars%nmo,
$ g_work,1,nbf_ao,1,pars%nmo)
call sym_movecs_adapt(ao_bas_han, 1.0d-8, g_work,
$ int_mb(pars%k_irs(i)),nmix(i))
if(.not.ga_destroy(g_work)) call errquit
$ ('bse_init: failed to destroy g_work',0, GA_ERR)
enddo
!Clean previous grid
call grid_cleanup(.false.)
!Print expected memory requirements
!Compute three-center integrals
call bse_ri_init(pars)
end subroutine bse_init

169
src/bsemol/bse_input.F Normal file
View file

@ -0,0 +1,169 @@
subroutine bse_input(rtdb)
c
c Input parser for molecular BSE module
c
c Daniel Mejia-Rodriguez
c November 2022
c
implicit none
c
#include "errquit.fh"
#include "util.fh"
#include "inp.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "stdio.fh"
c
integer rtdb
character*255 test
character(*), parameter :: pname = "bse_input: "
c
integer, parameter :: num_dirs = 13
character*11 dirs(num_dirs)
data dirs / 'method','tda','ncap','ncapr','singlet','triplet',
& 'outwin', 'nroots', 'maxroots', 'skipgw', 'maxiter',
& 'nstart', 'end' /
c
integer ind, itemp
double precision threshold, eta, temp
double precision, parameter :: ha2ev = 27.211386245988d0
c
c Set default values and put them in RTDB
c
call bse_defaults(rtdb)
c
c Parse the input file
c
10 if (.not. inp_read()) then
call errquit(pname//'failed to read input',0,INPUT_ERR)
endif
c
if (.not. inp_a(test)) then
call errquit(pname//'failed to read keyword',0,INPUT_ERR)
endif
c
if (.not.inp_match(num_dirs, .false., test, dirs, ind)) then
call errquit(pname//'unknown directive', 0, INPUT_ERR)
end if
c
goto (100, 200, 300, 301, 400, 500, 600, 700, 800, 900, 1000,
& 1100, 9999) ind
call errquit(pname//'unknown directive',ind, INPUT_ERR)
c
c -----------------------------------------------------------------
c
c method
c
100 if (.not. inp_a(test))
& call errquit(pname//'could not read method', 100, INPUT_ERR)
if (inp_compare(.false.,'analytic',test)) then
if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 101, RTDB_ERR)
elseif (inp_compare(.false.,'davidson',test)) then
if (.not.rtdb_put(rtdb, 'bse:davidson', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 102, RTDB_ERR)
if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .false.))
& call errquit(pname//'rtdb_put failed', 101, RTDB_ERR)
elseif (inp_compare(.false.,'lanczos',test)) then
if (.not.rtdb_put(rtdb, 'bse:lanczos', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 103, RTDB_ERR)
if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .false.))
& call errquit(pname//'rtdb_put failed', 101, RTDB_ERR)
else
call errquit(pname//'unknown BSE method',104,INPUT_ERR)
endif
goto 10
c
c use Tamm-Dancoff approximation
c
200 if (.not.rtdb_put(rtdb, 'bse:tda', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 200, RTDB_ERR)
goto 10
c
c use NCAP DD correction for eigenvalues
c
300 if (.not.rtdb_put(rtdb, 'bse:ncap', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)
goto 10
c
c use NCAPR DD correction for eigenvalues
c
301 if (.not.rtdb_put(rtdb, 'bse:ncapr', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 301, RTDB_ERR)
goto 10
c
c compute singlet transitions
c
400 if (.not.rtdb_put(rtdb, 'bse:singlet', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)
if (.not.rtdb_put(rtdb, 'bse:triplet', mt_log, 1, .false.))
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)
goto 10
c
c compute triplet transitions
c
500 if (.not.rtdb_put(rtdb, 'bse:singlet', mt_log, 1, .false.))
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)
if (.not.rtdb_put(rtdb, 'bse:triplet', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)
goto 10
c
c output energy window
c
600 if (.not.rtdb_put(rtdb, 'bse:dowindow', mt_log, 1, .true.))
& call errquit(pname//'rtdb_put failed', 600, RTDB_ERR)
if(.not.inp_f(temp))
& call errquit(pname//'could not read lower energy',600,INPUT_ERR)
if(.not.rtdb_put(rtdb,'bse:elower',mt_dbl,1,temp))
& call errquit(pname//'rtdb_put failed', 600, RTDB_ERR)
if(.not.inp_f(temp))
& call errquit(pname//'could not read upper energy',600,INPUT_ERR)
if(.not.rtdb_put(rtdb,'bse:eupper',mt_dbl,1,temp))
& call errquit(pname//'rtdb_put failed', 600, RTDB_ERR)
goto 10
c
c number of roots
c
700 if(.not.inp_i(itemp))
& call errquit(pname//'could not read no. of roots',700,INPUT_ERR)
if(.not.rtdb_put(rtdb,'bse:nroots',mt_int,1,itemp))
& call errquit(pname//'rtdb_put failed', 700, RTDB_ERR)
goto 10
c
c maximum number of trial vectors
c
800 if(.not.inp_i(itemp))
& call errquit(pname//'could not read maxroots',800,INPUT_ERR)
if(.not.rtdb_put(rtdb,'bse:nmax',mt_int,1,itemp))
& call errquit(pname//'rtdb_put failed', 800, RTDB_ERR)
goto 10
c
c Skip GW calculation
c
900 if(.not.rtdb_put(rtdb,'bse:skipgw',mt_log,1,.true.))
& call errquit(pname//'rtdb_put failed', 900, RTDB_ERR)
goto 10
c
c Maximum number of iterations
c
1000 if(.not.inp_i(itemp))
& call errquit(pname//'could not read maxiter', 1000, INPUT_ERR)
if(.not.rtdb_put(rtdb,'bse:maxiter',mt_int,1,itemp))
& call errquit(pname//'rtdb_put failed', 1000, RTDB_ERR)
goto 10
c
c Size of starting space
c
1100 if(.not.inp_i(itemp))
& call errquit(pname//'could not read nstart', 1100, INPUT_ERR)
if(.not.rtdb_put(rtdb,'bse:nspace', mt_int, 1, itemp))
& call errquit(pname//'rtdb_put failed', 1100, RTDB_ERR)
goto 10
c
c normal termination
c
9999 return
c
end subroutine bse_input

9
src/bsemol/bse_lanczos.F Normal file
View file

@ -0,0 +1,9 @@
subroutine bse_lanczos(pars)
implicit none
#include "bse.fh"
type(bse_params_t) :: pars
end subroutine bse_lanczos

515
src/bsemol/bse_oscstr.F Normal file
View file

@ -0,0 +1,515 @@
subroutine bse_oscstr(pars,omega,totpoles)
implicit none
#include "errquit.fh"
#include "stdio.fh"
#include "global.fh"
#include "mafdecls.fh"
#include "bse.fh"
#include "cdft.fh"
#include "rtdb.fh"
#include "sym.fh"
#include "util.fh"
#include "util_params.fh"
#include "inp.fh"
type(bse_params_t) :: pars
character(*), parameter :: pname = 'bse_oscstr: '
character(len=7) :: key
integer totpoles,npoles
double precision omega(*)
double precision oscstr(3),oscstr_tot,tmp,moms(20),maxr,r(1)
integer llpole,ulpole,iroot,ixyz,ipole,mag_transdip(3)
integer isym,jsym,i,j,k,l,pol,occ,vir,ioff
integer g_td(2),g_xy(2),g_dipmag, g_ovlp, g_work, g_corr
integer g_tdtot(2)
integer k_corr, l_corr, ulloop
character*4 oname,vname
integer g_dipx, g_dipy, g_dipz
integer g_qdrxx, g_qdrxy, g_qdrxz, g_qdryy, g_qdryz, g_qdrzz
integer g_work2, g_work3, g_x(2), g_y(2)
logical ltransden,stat
double precision origin(3)
double precision na,nb,s2,s2_full
double precision espread
data origin/0d0,0d0,0d0/
character*(nw_max_path_len) fn_civecs, fn_transden
logical,external :: dmat_to_file
integer,external :: ga_create_atom_blocked
character*5 spin(2)
double precision x(totpoles),y(totpoles)
double precision thresh ! Threshold for printing
parameter (thresh=5.0d-2)
double precision xe(3), xh(3), xe2_exp, xh2_exp, xe_exp_2
double precision xh_exp_2, xhxe_exp, dhtoe, dexc, sigma_h, sigma_e
double precision omega_norm, cov, pcc, dcd
spin(1)='alpha'
spin(2)='beta '
if(.not.rtdb_get(pars%rtdb,'bse:ltransden',mt_log,1,ltransden))
$ ltransden = .false.
if(.not.nga_create(mt_dbl,3,(/nbf_ao,nbf_ao,3/),'mag_dipole',
$ (/nbf_ao,-1,-1/),g_dipmag))
$ call errquit('bse_oscstr: nga_create failed g_dipmag',0,GA_ERR)
call ga_zero(g_dipmag)
call int_giao_1ega(ao_bas_han,ao_bas_han,g_dipmag,'angmom',origin,
$ 1,.false.)
call ga_scale(g_dipmag,0.5d0)
g_ovlp = ga_create_atom_blocked(geom, ao_bas_han, 'Overlap')
call ga_zero(g_ovlp)
call int_1e_ga(ao_bas_han, ao_bas_han, g_ovlp, 'overlap', oskel)
if(oskel) call sym_symmetrize(geom, ao_bas_han, .false., g_ovlp)
if (pars%ipol.eq.2) then
if(.not.ma_push_get(mt_dbl,pars%nmo*pars%nmo,'corr',
$ l_corr,k_corr))
$ call errquit('bse_oscstr: failed to create MA',0,MA_ERR)
if(.not.ga_create(mt_dbl,nbf_ao,pars%nmo,'work',0,0,g_work))
$ call errquit('bse_oscstr: failed to create GA',0,GA_ERR)
if(.not.ga_create(mt_dbl,pars%nmo,pars%nmo,'work',
$ pars%nmo,0,g_corr))
$ call errquit('bse_oscstr: failed to create GA',0,GA_ERR)
call ga_dgemm('n','n',nbf_ao,pars%nmo,nbf_ao,1d0,g_ovlp,
$ pars%g_movecs(2),0d0,g_work)
call ga_dgemm('t','n',pars%nmo,pars%nmo,nbf_ao,1d0,
$ pars%g_movecs(1),g_work,0d0,g_corr)
call ga_get(g_corr,1,pars%nmo,1,pars%nmo,dbl_mb(k_corr),
$ pars%nmo)
stat = ga_destroy(g_work).and.ga_destroy(g_corr)
if(.not.stat)
$ call errquit('bse_osctr: failed to destroy GA',58,GA_ERR)
endif
g_dipx = ga_create_atom_blocked(geom, ao_bas_han, 'X')
g_dipy = ga_create_atom_blocked(geom, ao_bas_han, 'Y')
g_dipz = ga_create_atom_blocked(geom, ao_bas_han, 'Z')
g_qdrxx = ga_create_atom_blocked(geom, ao_bas_han, 'XX')
g_qdrxy = ga_create_atom_blocked(geom, ao_bas_han, 'XY')
g_qdrxz = ga_create_atom_blocked(geom, ao_bas_han, 'XZ')
g_qdryy = ga_create_atom_blocked(geom, ao_bas_han, 'YY')
g_qdryz = ga_create_atom_blocked(geom, ao_bas_han, 'YZ')
g_qdrzz = ga_create_atom_blocked(geom, ao_bas_han, 'ZZ')
g_work = ga_create_atom_blocked(geom, ao_bas_han, 'work1')
g_work2 = ga_create_atom_blocked(geom, ao_bas_han, 'work2')
g_work3 = ga_create_atom_blocked(geom, ao_bas_han, 'work3')
call int_dip_ga(ao_bas_han, ao_bas_han, g_dipx, g_dipy, g_dipz)
call int_qdr_ga(ao_bas_han, ao_bas_han, g_qdrxx, g_qdrxy, g_qdrxz,
$ g_qdryy, g_qdryz, g_qdrzz)
npoles = pars%npoles(1)
if (pars%ipol.gt.1) npoles = npoles + pars%npoles(2)
if (pars%singlet) then
key = 'singlet'
else
key = 'triplet'
endif
! Get the lower and upper limits
llpole = -1
if (pars%dowindow .and. (.not.pars%davidson)) then
iroot = 0
do ipole=1,npoles
if (omega(ipole).lt.pars%elower) cycle
if (llpole.eq.-1) llpole = ipole
iroot = iroot + 1
if (omega(ipole).gt.pars%eupper) then
ulpole = ipole - 1
exit
endif
if (iroot.eq.pars%nroots) then
ulpole = ipole
exit
endif
enddo
else
llpole = 1
ulpole = pars%nroots
endif
! Get symmetry of ground state
isym = 1
if (pars%ipol.eq.2) then
do i=1,pars%ipol
do j=1,pars%nocc(i)
isym = ieor(isym-1,int_mb(pars%k_irs(i)+j-1)-1)+1
enddo
enddo
endif
call sym_irrepname(geom, isym, oname)
! Allocate transition density matrices
do i=1,pars%ipol
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'Transition density',
$ -1,-1,g_td(i))) call errquit
$ ('bse_oscstr: failed to create g_td',0,GA_ERR)
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'Transition density',
$ -1,-1,g_tdtot(i))) call errquit
$ ('bse_oscstr: failed to create g_tdtot',0,GA_ERR)
if(.not.ga_create(mt_dbl,pars%npoles(i),1,'X+Y vector',-1,1,
$ g_xy(i))) call errquit
$ ('bse_oscstr: failed to create g_xy',0,GA_ERR)
if(.not.ga_create(mt_dbl,pars%npoles(i),1,'X vector',-1,1,
$ g_x(i))) call errquit
$ ('bse_oscstr: failed to create g_x',0,GA_ERR)
if(.not.ga_create(mt_dbl,pars%npoles(i),1,'Y vector',-1,1,
$ g_y(i))) call errquit
$ ('bse_oscstr: failed to create g_y',0,GA_ERR)
enddo
! Compute the oscillator strengths around an energy window
iroot = 0
do ipole=llpole,ulpole
iroot = iroot + 1
oscstr_tot = 0d0
moms(:) = 0d0
oscstr(:) = 0d0
mag_transdip = 0d0
call ga_copy_patch('n',pars%g_apb,1,pars%npoles(1),ipole,ipole,
& g_x(1),1,pars%npoles(1),1,1)
if(pars%ipol.gt.1) then
call ga_copy_patch('n',pars%g_apb,pars%npoles(1)+1,npoles,
& ipole,ipole,g_x(2),1,pars%npoles(2),1,1)
endif
call tddft_transfm(1, g_x, pars%g_movecs, nbf_ao, pars%nocc,
& (/pars%nmo, pars%nmo/), (/0,0/), (/0,0/),
& pars%ipol, g_td, 1)
call ga_transpose(g_td(1), g_tdtot(1))
if(pars%ipol.gt.1) call ga_transpose(g_td(2), g_tdtot(2))
if(pars%davidson) then
call ga_copy_patch('n',pars%g_amb,1,pars%npoles(1),
I ipole,ipole,
& g_y(1),1,pars%npoles(1),1,1)
if(pars%ipol.gt.1) then
call ga_copy_patch('n',pars%g_amb, pars%npoles(1)+1, npoles,
& ipole, ipole, g_y(2), 1, pars%npoles(2),
& 1, 1)
endif
call tddft_transfm(1, g_y, pars%g_movecs, nbf_ao, pars%nocc,
& (/pars%nmo, pars%nmo/), (/0,0/), (/0,0/),
& pars%ipol, g_td, 1)
call ga_add(1.0d0, g_td(1), 1.0d0, g_tdtot(1), g_td(1))
if(pars%ipol.gt.1)
& call ga_add(1.0d0, g_td(2), 1.0d0, g_tdtot(2), g_td(2))
else
call ga_copy(g_tdtot(1), g_td(1))
if (pars%ipol.gt.1) call ga_copy(g_tdtot(2), g_td(2))
endif
do i=1,pars%ipol
call multipole_density(ao_bas_han,origin,3,g_td(i),moms,20)
do ixyz=1,3
mag_transdip(ixyz) = mag_transdip(ixyz) +
$ nga_ddot_patch(g_td(i),'n',(/1,1/),(/nbf_ao,nbf_ao/),
$ g_dipmag,'n',(/1,1,ixyz/),
$ (/nbf_ao,nbf_ao,ixyz/))
enddo
enddo
if(pars%ipol.gt.1) call ga_add(1d0,g_td(1),1d0,g_td(2),g_td(1))
! Indicators
call ga_dgemm('t','n',nbf_ao,nbf_ao,nbf_ao,1d0,g_td(1),g_ovlp,
$ 0d0,g_work)
call ga_dgemm('n','n',nbf_ao,nbf_ao,nbf_ao,1d0,g_work,g_td(1),
$ 0d0,g_work2)
omega_norm = ga_ddot(g_work2,g_ovlp)
xe(1) = ga_ddot(g_work2,g_dipx)/omega_norm
xe(2) = ga_ddot(g_work2,g_dipy)/omega_norm
xe(3) = ga_ddot(g_work2,g_dipz)/omega_norm
xe_exp_2 = xe(1)**2 + xe(2)**2 + xe(3)**2
xe2_exp = (ga_ddot(g_work2,g_qdrxx) +
$ ga_ddot(g_work2,g_qdryy) +
$ ga_ddot(g_work2,g_qdrzz))/omega_norm
call ga_dgemm('n','n',nbf_ao,nbf_ao,nbf_ao,1d0,g_td(1),g_ovlp,
$ 0d0,g_work)
call ga_dgemm('n','t',nbf_ao,nbf_ao,nbf_ao,1d0,g_work,g_td(1),
$ 0d0,g_work2)
xh(1) = ga_ddot(g_work2,g_dipx)/omega_norm
xh(2) = ga_ddot(g_work2,g_dipy)/omega_norm
xh(3) = ga_ddot(g_work2,g_dipz)/omega_norm
xh_exp_2 = xh(1)**2 + xh(2)**2 + xh(3)**2
xh2_exp = (ga_ddot(g_work2,g_qdrxx) +
$ ga_ddot(g_work2,g_qdryy) +
$ ga_ddot(g_work2,g_qdrzz))/omega_norm
xhxe_exp = 0d0
call ga_dgemm('n','n',nbf_ao,nbf_ao,nbf_ao,1d0,g_td(1),g_dipx,
$ 0d0,g_work)
call ga_dgemm('n','t',nbf_ao,nbf_ao,nbf_ao,1d0,g_work,g_td(1),
$ 0d0,g_work2)
xhxe_exp = xhxe_exp + ga_ddot(g_work2,g_dipx)/omega_norm
call ga_dgemm('n','n',nbf_ao,nbf_ao,nbf_ao,1d0,g_td(1),g_dipy,
$ 0d0,g_work)
call ga_dgemm('n','t',nbf_ao,nbf_ao,nbf_ao,1d0,g_work,g_td(1),
$ 0d0,g_work2)
xhxe_exp = xhxe_exp + ga_ddot(g_work2,g_dipy)/omega_norm
call ga_dgemm('n','n',nbf_ao,nbf_ao,nbf_ao,1d0,g_td(1),g_dipz,
$ 0d0,g_work)
call ga_dgemm('n','t',nbf_ao,nbf_ao,nbf_ao,1d0,g_work,g_td(1),
$ 0d0,g_work2)
xhxe_exp = xhxe_exp + ga_ddot(g_work2,g_dipz)/omega_norm
sigma_h = dsqrt(xh2_exp - xh_exp_2)
sigma_e = dsqrt(xe2_exp - xe_exp_2)
dhtoe = dsqrt((xe(1)-xh(1))**2 +
$ (xe(2)-xh(2))**2 +
$ (xe(3)-xh(3))**2)
dexc = dsqrt(xe2_exp + xh2_exp - 2d0*xhxe_exp)
cov = xhxe_exp - xe(1)*xh(1) - xe(2)*xh(2) - xe(3)*xh(3)
pcc = cov/sigma_e/sigma_h
dcd = dhtoe - 0.5d0*(sigma_e + sigma_h)
if (ltransden) then
call util_file_name('tdens',.false.,.false.,fn_transden)
call createuniquefilename(fn_transden,ipole)
call util_file_name_resolve(fn_transden,.false.)
if (pars%ipol.gt.1)
& call ga_add(1d0,g_td(1),1d0,g_td(2),g_td(1))
call ga_symmetrize(g_td(1))
if(.not. dmat_to_file(g_td(1),fn_transden))
& call errquit('bse_oscstr: dmat_to_file error',0,CALC_ERR)
endif
if (pars%ipol.eq.1) then
do i=1,20
moms(i) = moms(i)*dsqrt(2.0d0)
enddo
endif
if (pars%davidson) then
call ga_get(pars%g_apb,1,npoles,ipole,ipole,x,npoles)
call ga_get(pars%g_amb,1,npoles,ipole,ipole,y,npoles)
endif
if (pars%davidson.and.pars%ipol.eq.2) then
call tddft_s2(pars%tda, (/0,0/), pars%nocc,
& (/pars%nmo,pars%nmo/), (/0,0/),
& pars%npoles,x,x(pars%npoles(1)+1),y,
& y(pars%npoles(1)+1),dbl_mb(k_corr),s2)
call bse_s2(pars%nocc,pars%nmo,pars%npoles,x,
& x(pars%npoles(1)+1), y, y(pars%npoles(1)+1), dbl_mb(k_corr),
& s2_full)
if (abs(s2_full).lt.0.25d0) then
key = 'singlet'
elseif (abs(abs(s2_full)-2.0d0).lt.0.25) then
key = 'triplet'
else
key = ' '
endif
endif
if (pars%me.eq.0) then
call tddft_oscstr(moms, omega(ipole), mag_transdip,
& oscstr, oscstr_tot)
l = 0
maxr = 0.0d0
pol = 1
occ = 1
vir = pars%nocc(1) + 1
do i=1,pars%ipol
if (pars%npoles(i).eq.0) cycle
do j=1,pars%nocc(i)
do k=pars%nocc(i)+1,pars%nmo
l = l + 1
call ga_get(pars%g_apb,l,l,ipole,ipole,r,1)
if (dabs(r(1)).gt.maxr) then
maxr = dabs(r(1))
pol = i
occ = j
vir = k
endif
enddo
enddo
enddo
jsym = ieor(isym-1,int_mb(pars%k_irs(pol)+occ-1)-1)+1
jsym = ieor(jsym-1,int_mb(pars%k_irs(pol)+vir-1)-1)+1
call sym_irrepname(geom, jsym, vname)
write(luout,9110) iroot,key,vname,omega(ipole),
& omega(ipole)*ha2ev
if (pars%davidson.and.pars%ipol.eq.2) then
write(luout,*)
write(luout,9180) s2
write(luout,9190) s2_full
endif
if (pars%singlet) then
write(luout,9200) moms(2),moms(3),moms(4),
& moms(5),moms(6),moms(7),
& moms(8),moms(9),moms(10),
& oscstr(1),oscstr(2),oscstr(3),
& oscstr_tot
else
write(luout,9210)
endif
write(luout,*)
write(luout,9340) omega_norm
write(luout,9300) sigma_e*cau2ang
write(luout,9310) sigma_h*cau2ang
write(luout,9320) dhtoe*cau2ang
write(luout,9330) dexc*cau2ang
write(luout,9360) dcd*cau2ang
write(luout,9350) pcc
if (pars%davidson) then
write(luout,*)
if (ipol.eq.1) then
l=0
do j=1,pars%nocc(1)
do k=pars%nocc(1)+1,pars%nmo
l=l+1
if (pars%tda) then
r=x(l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9120) j,oname,k,vname,r
else
r=x(l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9130) j,oname,k,vname,r
r=y(l)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9140) j,oname,k,vname,r
endif
enddo
enddo
else
do i=1,ipol
ioff = (i-1)*pars%npoles(1)
l=0
do j=1,pars%nocc(i)
do k=pars%nocc(i)+1,pars%nmo
l=l+1
if (pars%tda) then
r = x(ioff+l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9150) j,spin(i),oname,
3 k,spin(i),vname,r
else
r = x(ioff+l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9160) j,spin(i),oname,
3 k,spin(i),vname,r
r=y(ioff+l)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9170) j,spin(i),oname,
3 k,spin(i),vname,r
endif
enddo
enddo
enddo
endif
call util_flush(LuOut)
endif
endif
call ga_sync()
enddo
if (.not.ga_destroy(g_dipmag))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
do i=1,pars%ipol
if (.not.ga_destroy(g_td(i)))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
if (.not.ga_destroy(g_xy(i)))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
if (.not.ga_destroy(g_tdtot(i)))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
if (.not.ga_destroy(g_x(i)))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
if (.not.ga_destroy(g_y(i)))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
enddo
if(.not.(ga_destroy(g_dipx).and.ga_destroy(g_dipy).and.
$ ga_destroy(g_dipz).and.ga_destroy(g_qdrxx).and.
$ ga_destroy(g_qdrxy).and.ga_destroy(g_qdrxz).and.
$ ga_destroy(g_qdryy).and.ga_destroy(g_qdryz).and.
$ ga_destroy(g_qdrzz).and.ga_destroy(g_ovlp).and.
$ ga_destroy(g_work).and.ga_destroy(g_work2).and.
$ ga_destroy(g_work3)))
& call errquit('bse_oscstr: ga_destroy failed',0,GA_ERR)
9110 format(2x,
1 '---------------------------------------------------------------
3-------------',
2 /,2x,'Root',i4,1x,a7,1x,a4,f22.9,' a.u.',f22.4,' eV ',/,2x,
1 '---------------------------------------------------------------
3-------------')
9200 format(5x,'Transition Moments X',f9.5,' Y',f9.5,' Z',f9.5
1 ,/,5x,'Transition Moments XX',f9.5,' XY',f9.5,' XZ',f9.5
2 ,/,5x,'Transition Moments YY',f9.5,' YZ',f9.5,' ZZ',f9.5
5 ,/,5x,'Dipole Oscillator Strength',18x,f14.10
5 ,/,5x,'Electric Quadrupole ',18x,f14.10
5 ,/,5x,'Magnetic Dipole ',18x,f14.10
5 ,/,5x,'Total Oscillator Strength ',18x,f14.10)
9210 format(5x,'Transition Moments Spin forbidden'
1 ,/,5x,'Oscillator Strength Spin forbidden')
9120 format(5x,'Occ.',i5,2x,a4,'--- Virt.',i5,2x,a4,f12.5)
9130 format(5x,'Occ.',i5,2x,a4,'--- Virt.',i5,2x,a4,f8.5,' X')
9140 format(5x,'Occ.',i5,2x,a4,'--- Virt.',i5,2x,a4,f8.5,' Y')
9150 format(5x,'Occ.',i5,1x,a5,1x,a4,'--- Virt.',
1 i5,1x,a5,1x,a4,f12.5)
9160 format(5x,'Occ.',i5,1x,a5,1x,a4,'--- Virt.',
1 i5,1x,a5,1x,a4,f8.5,' X')
9170 format(5x,'Occ.',i5,1x,a5,1x,a4,'--- Virt.',
1 i5,1x,a5,1x,a4,f8.5,' Y')
9180 format(2x,' <S2> = ',f8.4)
9190 format(2x,' D<S2> = ',f8.4)
9300 format(5x,' sigma_e = ',f8.4, ' Angstrom')
9310 format(5x,' sigma_h = ',f8.4, ' Angstrom')
9320 format(5x,' |d_h->e| = ',f8.4, ' Angstrom')
9330 format(5x,' d_exc = ',f8.4, ' Angstrom')
9360 format(5x,' d_CD = ',f8.4, ' Angstrom')
9350 format(5x,'Pearson C.C. = ',f8.4)
9340 format(5x,'Exciton norm = ',f8.4)
end

90
src/bsemol/bse_puteris.F Normal file
View file

@ -0,0 +1,90 @@
subroutine bse_puteris(erimo,nmo,nstates,nri,nri_me,rilo,ilower,
$ iupper,jlower,jupper,eri,blk,polelo,
$ polehi,mynpoles)
implicit none
#include "mafdecls.fh"
#include "errquit.fh"
#include "global.fh"
character(len=2),intent(in) :: blk
integer,intent(in) :: nmo, nstates, nri, nri_me, rilo
integer,intent(in) :: ilower, iupper, jlower, jupper
integer,intent(out) :: eri,polelo,polehi,mynpoles
double precision,intent(in),dimension(nstates,nmo,nri_me) :: erimo
character(len=12),parameter :: pname = 'gw_puteris: '
character(len=6) :: erilabel,idxlabel
integer ilo,ihi,jlo,jhi,imo,jmo,ipair,npairs,me,rihi
integer ni,nj,l_tmp,k_tmp,k_loc,ld
integer l_pair,k_pair,mynpairs,jpair,iproc,nprocs,ulpair,llpair
integer ibatch,nbatch,myllpair,myulpair,ma_dbl
integer lmap, kmap
logical ok
integer,external :: gw_get_pair
integer l_oo,k_oo,l_ov,k_ov,l_vv,k_vv
double precision temp(nri)
ni = iupper - ilower + 1
if (ni.lt.1) return
nj = jupper - jlower + 1
if (nj.lt.1) return
rihi = rilo + nri_me - 1
me = ga_nodeid()
nprocs = ga_nnodes()
erilabel = 'eri '//blk
idxlabel = 'idx '//blk
npairs = ni*nj
mynpairs = npairs/nprocs + 1
if (.not.ga_create(mt_dbl,nri,npairs,erilabel,nri,0,eri))
$ call errquit(pname//'can''t create '//blk//' block',0,GA_ERR)
ma_dbl = ma_inquire_stack(mt_dbl)/nri_me
nbatch = mynpairs/ma_dbl + 1
mynpairs = npairs/nbatch
if (mod(npairs,nbatch).ne.0) mynpairs = mynpairs + 1
mynpairs = min(npairs,mynpairs)
if (.not.ma_push_get(mt_dbl,mynpairs*nri_me,'temp',l_pair,k_pair))
$ call errquit(pname//'can''t create temporary array',0,MA_ERR)
ulpair = 0
do iproc=0,nprocs-1
llpair = ulpair + 1
ulpair = ulpair + npairs/nprocs
if (iproc.lt.mod(npairs,nprocs)) ulpair = ulpair + 1
if (llpair.gt.ulpair) goto 101
do ibatch=1,nbatch
myllpair = (ibatch-1)*mynpairs + llpair
myulpair = min(myllpair + mynpairs - 1,ulpair)
do ipair=myllpair,myulpair
jmo = (ipair-1)/ni + 1
imo = ipair - (jmo-1)*ni + ilower - 1
jmo = jmo + jlower - 1
call dcopy(nri_me,erimo(imo,jmo,1),nstates*nmo,
$ dbl_mb(k_pair+(ipair-myllpair)*nri_me),1)
enddo
call ga_put(eri,rilo,rihi,myllpair,myulpair,dbl_mb(k_pair),
$ nri_me)
enddo
101 continue
call ga_sync()
enddo
ok = ma_chop_stack(l_pair)
call ga_distribution(eri,me,ld,ld,polelo,polehi)
if (blk .eq. 'ov') then
mynpoles = polehi - polelo + 1
endif
end subroutine

327
src/bsemol/bse_ri_init.F Normal file
View file

@ -0,0 +1,327 @@
subroutine bse_ri_init(pars)
c
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "cdft.fh"
#include "bas.fh"
#include "bse.fh"
#include "stdio.fh"
#include "util.fh"
#include "global.fh"
#ifdef SCALAPACK
#include "scaleig.fh"
#endif
c
type(bse_params_t) :: pars
c
character(len=*), parameter :: pname = "bse_ri_init: "
c
logical iolgc
integer g_2ceris, g_2cinv, info, lmiss
c
integer ilo,ihi,jlo,jhi
integer lbuf,lscr,ltmp,leri,lvecs
integer l_buf,k_buf,l_scr,k_scr,l_tmp,k_tmp,k_eri,k_mo,l_mo
integer l_vecs,k_vecs,l_idx,k_idx
integer l_shlbuf, k_shlbuf, l_mobuf, k_mobuf
integer me, isp, addr
integer g_3ceris,shlo,shhi,offset,ld,g_tmp
integer nmo,nvir,rihi,rilo,nri,nri_me,nocc,offk
integer ish,jsh,ksh,ni,nj,nk,nij
integer ifirst,jfirst,kfirst
integer ilast,jlast,klast,imaxthreads
integer aopairs,mopairs(2),idum,totmopairs
logical IeqJ,ok,oactive
double precision eri_est,timer,timing
c
double precision schwarz_shell
external schwarz_shell
c
#ifdef USE_OPENMP
integer,external :: omp_get_max_threads
#endif
timer = util_wallsec()
c
c Create Global Arrays
c
if (.not.ga_create(mt_dbl,nbf_cd,nbf_cd,'2c ERI',
$ 0, nbf_cd, g_2ceris))
$ call errquit(pname//'Error creating 2c ERI',0,GA_ERR)
if (.not.ga_create(mt_dbl,nbf_cd,nbf_cd,'2c INV',
$ 0, nbf_cd, g_2cinv))
$ call errquit(pname//'Error creating 2c ERI',0,GA_ERR)
call ga_zero(g_2ceris)
c
c Compute 2-center ERIs
c
if (pars%me.eq.0) write(luout,9000)
9000 format(/,10x,'Computing 2-center integrals')
c
call dft_get2eri(cd_bas_han, g_2ceris, .false.)
c
c Obtain inverse Cholesky factor
c
if (pars%me.eq.0) write(luout,9010)
9010 format(10x,'Computing Inverse Cholesky factor')
c
call ga_sync()
call ga_chol_inv_seq(g_2ceris, g_2cinv, info)
if (info.ne.0) then
if (pars%me.eq.0) then
write(luout,*)
write(luout,*)' Problem in performing Cholesky'
write(luout,*)' Obtaining square root via truncated '
write(luout,*)' eigenvalue decomposition '
write(luout,*)
endif
call dft_invdiag(g_2ceris, g_2cinv, nbf_cd)
endif
c
if (.not. ga_destroy(g_2ceris))
$ call errquit(pname//'failed to destrou g_2ceri',0,GA_ERR)
c
c Allocate storage
c
nri = nbf_cd
nmo = pars%nmo
aopairs = (nbf_ao*(nbf_ao+1))/2
do isp=1,pars%ipol
mopairs(isp) = pars%nmo*pars%nmo
enddo
totmopairs = sum(mopairs(1:pars%ipol))
if (.not.ga_create(mt_dbl,1,nri,'3c ERI',aopairs,0,
$ g_3ceris))
$ call errquit(pname//'can''t create 3c ERIs',0,GA_ERR)
call ga_distribution(g_3ceris,pars%me,ilo,ihi,rilo,rihi)
nri_me = rihi - rilo + 1
ok = ga_destroy(g_3ceris)
c call ga_access(g_3ceris,ilo,ihi,rilo,rihi,k_eri,ld)
c call dfill(aopairs*nri_me,0d0,dbl_mb(k_eri),1)
shlo = 0
do ksh=1,nshells_cd
ok = bas_cn2bfr(cd_bas_han,ksh,kfirst,klast)
if (klast.lt.rilo) cycle
if (shlo.eq.0) then
offset = (rilo-kfirst)
shlo = ksh
endif
shhi = ksh
if (klast.ge.rihi) exit
enddo
oactive = nri_me.gt.0
call int_mem_2e3c(lbuf, lscr)
if (.not.ma_alloc_get(mt_dbl,lbuf, 'RI buffer',l_buf, k_buf))
$ call errquit(pname//'can''t get buffer space',0,MA_ERR)
lscr = max(lscr,nbf_cd*nbf_ao_mxnbf_cn**2)
if (.not.ma_alloc_get(mt_dbl,lscr, 'RI scratch', l_scr, k_scr))
$ call errquit(pname//'can''t get scratch space',0,MA_ERR)
if (.not.ma_alloc_get(mt_dbl,nri_me*nbf_ao_mxnbf_cn**2,
$ 'Shell buffer',l_shlbuf, k_shlbuf))
$ call errquit(pname//'can''t allocate shell buffer',0,MA_ERR)
if (.not.ma_alloc_get(mt_dbl,
$ nri_me*nbf_ao_mxnbf_cn*nmo*pars%ipol,
$ 'MO buffer',l_mobuf, k_mobuf))
$ call errquit(pname//'can''t allocate shell buffer',0,MA_ERR)
c
c Compute three-center ERIs in MO representation
c
call ga_sync()
c
if (pars%me.eq.0) write(luout,9020)
9020 format(10x,'Computing 3-center integrals')
lvecs = nbf_ao*nmo*pars%ipol
if (.not.ma_alloc_get(mt_dbl,lvecs,'MOVECS',l_vecs,k_vecs))
$ call errquit(pname//'can''t get movecs space',0,MA_ERR)
if (.not.ma_alloc_get(mt_dbl,totmopairs*nri_me,'MO Eris',
$ l_mo,k_mo))
$ call errquit(pname//'can''t get MO eris space',0,MA_ERR)
call dfill(totmopairs*nri_me,0.0d0,dbl_mb(k_mo),1)
do isp=1,pars%ipol
call ga_get(pars%g_movecs(isp),1,nbf_ao,1,nmo,
$ dbl_mb(k_vecs+(isp-1)*nbf_ao*nmo),nbf_ao)
enddo
c
#ifdef USE_OPENMP
iMaxThreads = omp_get_max_threads()
call util_blas_set_num_threads(iMaxThreads)
#endif
do ish=1,nshells_ao
ok = bas_cn2bfr(ao_bas_han,ish,ifirst,ilast)
ni = ilast - ifirst + 1
call dfill(ni*nmo*pars%ipol*nri_me,0d0,dbl_mb(k_mobuf),1)
do jsh=1,nshells_ao
eri_est = schwarz_shell(ish,jsh)
if (eri_est.lt.pars%tol2e) cycle
IeqJ = ish.eq.jsh
ok = bas_cn2bfr(ao_bas_han,jsh,jfirst,jlast)
nj = jlast - jfirst + 1
nij = ni*nj
do ksh=shlo,shhi
ok = bas_cn2bfr(cd_bas_han, ksh, kfirst, klast)
call int_2e3c(cd_bas_han,ksh,ao_bas_han,ish,jsh,lscr,
$ dbl_mb(k_scr),lbuf,dbl_mb(k_buf))
if (ksh.eq.shlo) then
kfirst = rilo
offk = k_buf + offset*nij
else
offk = k_buf
endif
if (ksh.eq.shhi) klast = rihi
nk = klast - kfirst + 1
kfirst = kfirst - rilo + 1
call dcopy(ni*nj*nk,dbl_mb(offk),1,
$ dbl_mb(k_shlbuf+(kfirst-1)*ni*nj),1)
enddo
do isp=1,pars%ipol
call gw_mo2(dbl_mb(k_shlbuf),
$ dbl_mb(k_mobuf+(isp-1)*nmo*ni*nri_me),
$ dbl_mb(k_vecs+(isp-1)*nbf_ao*nmo),nri_me,
$ nbf_ao,ni,nj,nmo,nmo,ifirst,
$ jfirst,ieqj)
enddo
enddo
do isp=1,pars%ipol
call gw_mo(dbl_mb(k_mobuf+(isp-1)*nmo*ni*nri_me),
$ dbl_mb(k_mo+(isp-1)*mopairs(1)*nri_me),
$ dbl_mb(k_vecs+(isp-1)*nbf_ao*nmo),nri_me,
$ nbf_ao,ni,nmo,nmo,ifirst)
enddo
enddo
c
ok = ma_free_heap(l_buf) .and. ma_free_heap(l_scr) .and.
$ ma_free_heap(l_shlbuf) .and.
$ ma_free_heap(l_mobuf)
if (.not.ok)
$ call errquit(pname//'can''t free heap',165,MA_ERR)
c
c
#ifdef GWDEBUG
if (pars%me.eq.0) write(*,*) ' Transform ERIs'
#endif
call ga_sync()
ok = ma_free_heap(l_vecs)
if (.not.ok)
$ call errquit(pname//'can''t free heap',176,MA_ERR)
c
call ga_sync()
c
c Distribute and orthonormalize ERIs
c
timer = util_wallsec()
#ifdef GWDEBUG
if (pars%me.eq.0) write(*,*) ' Distribute ERIs'
#endif
addr = k_mo
do isp=1,pars%ipol
if (nmo.eq.0) cycle
if (isp.eq.2) addr = k_mo + mopairs(1)*nri_me
! OO block
call bse_puteris(dbl_mb(addr),nmo,nmo,nri,nri_me,
$ rilo,1,pars%nocc(isp),1,pars%nocc(isp),
$ pars%g_erioo(isp),'oo',pars%oolo(isp),
$ pars%oohi(isp),idum)
! OV block
call bse_puteris(dbl_mb(addr),nmo,nmo,nri,nri_me,
$ rilo,pars%nocc(isp)+1,nmo,1,pars%nocc(isp),
$ pars%g_eriov(isp),'ov',pars%ovlo(isp),
$ pars%ovhi(isp),pars%mynpoles(isp))
! VV block
call bse_puteris(dbl_mb(addr),nmo,nmo,nri,nri_me,
$ rilo,pars%nocc(isp)+1,nmo,
$ pars%nocc(isp)+1,nmo,pars%g_erivv(isp),
$ 'vv',pars%vvlo(isp),pars%vvhi(isp),idum)
enddo
ok = ma_free_heap(l_mo)
#ifdef GWDEBUG
if (pars%me.eq.0) write(*,*) ' Read inverse'
#endif
c
if (.not.ma_alloc_get(mt_dbl,nri**2,'RI temporary',l_tmp,k_tmp))
$ call errquit(pname//'can''t get temp space',0,MA_ERR)
if (pars%me.eq.0) then
call ga_get(g_2cinv,1,nbf_cd,1,nbf_cd,dbl_mb(k_tmp),nbf_cd)
endif
call ga_sync()
if (.not. ga_destroy(g_2cinv))
$ call errquit(pname//'failed to destroy g_2cinv',0,GA_ERR)
call ga_brdcst(1038,dbl_mb(k_tmp),ma_sizeof(mt_dbl,nri**2,
$ mt_byte),0)
#ifdef GWDEBUG
if (pars%me.eq.0) write(*,*) ' Orthogonalize ERIs'
#endif
do isp=1,pars%ipol
c
c Transform OO block
c
call ga_distribution(pars%g_erioo(isp),pars%me,ilo,ihi,jlo,jhi)
if (jlo.gt.jhi) goto 101
call ga_access(pars%g_erioo(isp),ilo,ihi,jlo,jhi,k_eri,ld)
nij = jhi-jlo+1
call dtrmm('l','l','n','n',nbf_cd,nij,1.0d0,dbl_mb(k_tmp),
$ nbf_cd,dbl_mb(k_eri),nbf_cd)
call ga_release_update(pars%g_erioo(isp),ilo,ihi,jlo,jhi)
c
c Transform OV block
c
101 continue
call ga_distribution(pars%g_eriov(isp),pars%me,ilo,ihi,jlo,jhi)
if (jlo.gt.jhi) goto 102
call ga_access(pars%g_eriov(isp),ilo,ihi,jlo,jhi,k_eri,ld)
nij = jhi-jlo+1
call dtrmm('l','l','n','n',nbf_cd,nij,1.0d0,dbl_mb(k_tmp),
$ nbf_cd,dbl_mb(k_eri),nbf_cd)
call ga_release_update(pars%g_eriov(isp),ilo,ihi,jlo,jhi)
c
c Transform VV block
c
102 continue
if (nmo.gt.pars%nocc(isp)) then
call ga_distribution(pars%g_erivv(isp),pars%me,ilo,ihi,
$ jlo,jhi)
if (jlo.gt.jhi) goto 103
call ga_access(pars%g_erivv(isp),ilo,ihi,jlo,jhi,k_eri,ld)
nij = jhi-jlo+1
call dtrmm('l','l','n','n',nbf_cd,nij,1.0d0,dbl_mb(k_tmp),
$ nbf_cd,dbl_mb(k_eri),nbf_cd)
call ga_release_update(pars%g_erivv(isp),ilo,ihi,jlo,jhi)
endif
103 continue
enddo
c
c Free heap
c
ok = ma_free_heap(l_tmp)
if (.not.ok)
$ call errquit(pname//'failed to free heap',224,MA_ERR)
c
#ifdef USE_OPENMP
call util_blas_set_num_threads(1)
#endif
return
end subroutine

99
src/bsemol/bse_s2.F Normal file
View file

@ -0,0 +1,99 @@
subroutine bse_s2(nocc, nmo, npoles, xup, xdw, yup, ydw,
$ ovlp, s2)
implicit none
integer nocc(2), npoles(2), nmo
double precision xup(npoles(1)), xdw(npoles(2))
double precision yup(npoles(1)), ydw(npoles(2))
double precision s2
double precision ovlp(nmo,nmo)
double precision ov,xx,yy
integer i,j,a,b,ipole,jpole,nvir(2)
integer llipole,ulipole,lljpole,uljpole
double precision, external :: ddot
s2 = 0d0
nvir(1) = nmo - nocc(1)
nvir(2) = nmo - nocc(2)
! Up-Up
do i=1,nocc(1)
do j=1,nocc(1)
llipole = (i-1)*nvir(1) + 1
ulipole = i*nvir(1)
lljpole = (j-1)*nvir(1) + 1
uljpole = j*nvir(1)
ov = dot_product(ovlp(i,1:nocc(2)),ovlp(j,1:nocc(2)))
xx = dot_product(xup(llipole:ulipole),xup(lljpole:uljpole))
yy = dot_product(yup(llipole:ulipole),yup(lljpole:uljpole))
s2 = s2 + ov*(xx+yy)
enddo
enddo
do a=nocc(1)+1,nmo
do b=nocc(1)+1,nmo
ov = dot_product(ovlp(a,1:nocc(2)),ovlp(b,1:nocc(2)))
llipole = a - nocc(1)
lljpole = b - nocc(1)
s2 = s2 - ov*ddot(nocc(1),xup(llipole),nvir(1),
$ xup(lljpole),nvir(1))
s2 = s2 - ov*ddot(nocc(1),yup(llipole),nvir(1),
$ yup(lljpole),nvir(1))
enddo
enddo
! Dw-Dw
do i=1,nocc(2)
do j=1,nocc(2)
llipole = (i-1)*nvir(2) + 1
ulipole = i*nvir(2)
lljpole = (j-1)*nvir(2) + 1
uljpole = j*nvir(2)
ov = dot_product(ovlp(1:nocc(1),i),ovlp(1:nocc(1),j))
xx = dot_product(xdw(llipole:ulipole),xdw(lljpole:uljpole))
yy = dot_product(ydw(llipole:ulipole),ydw(lljpole:uljpole))
s2 = s2 + ov*(xx+yy)
enddo
enddo
do a=nocc(2)+1,nmo
do b=nocc(2)+1,nmo
ov = dot_product(ovlp(1:nocc(1),a),ovlp(1:nocc(1),b))
llipole = a - nocc(2)
lljpole = b - nocc(2)
s2 = s2 - ov*ddot(nocc(2),xdw(llipole),nvir(2),
$ xdw(lljpole),nvir(2))
s2 = s2 - ov*ddot(nocc(2),ydw(llipole),nvir(2),
$ ydw(lljpole),nvir(2))
enddo
enddo
! Up-Dw
do ipole=1,npoles(1)
i = ipole/nvir(1) + 1
a = ipole - (i-1)*nvir(1) + nocc(1)
do jpole=1,npoles(2)
j = jpole/nvir(2) + 1
b = jpole - (j-1)*nvir(2) + nocc(2)
s2 = s2 -
$ 2d0*xup(ipole)*xdw(jpole)*ovlp(i,j)*ovlp(a,b) -
$ 2d0*yup(ipole)*ydw(jpole)*ovlp(i,j)*ovlp(a,b) +
$ 2d0*xup(ipole)*ydw(jpole)*ovlp(a,j)*ovlp(i,b)
enddo
enddo
! Dw-Up
do ipole=1,npoles(2)
i = ipole/nvir(2) + 1
a = ipole - (i-1)*nvir(2) + nocc(2)
do jpole=1,npoles(1)
j = jpole/nvir(1) + 1
b = jpole - (j-1)*nvir(1) + nocc(1)
s2 = s2 +
$ 2d0*xdw(ipole)*yup(jpole)*ovlp(j,a)*ovlp(b,i)
enddo
enddo
end subroutine

54
src/bsemol/bse_wmn.F Normal file
View file

@ -0,0 +1,54 @@
subroutine bse_wmn(pars,w,nmo,nri,ipol)
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "bse.fh"
type(bse_params_t) :: pars
integer nmo, nri, ipol
character(*), parameter :: pname = 'bse_wmn: '
double precision w(*)
integer isp,ilo,ihi,jlo,jhi,keri,ld,g_tmp
do isp=1,ipol
if (.not.pars%tda) then
! copy and transform occ-vir block
if(.not.ga_duplicate(pars%g_eriov(isp),pars%g_wov(isp),'Wov'))
& call errquit(pname//'could not duplicate ERI OV',0,GA_ERR)
call ga_copy(pars%g_eriov(isp),pars%g_wov(isp))
if(pars%mynpoles(isp).eq.0) cycle
call ga_access(pars%g_wov(isp),1,nri,pars%ovlo(isp),
& pars%ovhi(isp),keri,ld)
call dtfsm('n','l','l','n','n',nri,pars%mynpoles(isp),1d0,
& w,dbl_mb(keri),nri)
call ga_release_update(pars%g_wov(isp),1,nri,pars%ovlo(isp),
& pars%ovhi(isp))
endif
! transform occ-occ and vir-vir blocks in place
call ga_distribution(pars%g_erioo(isp),pars%me,ilo,ihi,jlo,jhi)
if (jlo.gt.jhi) goto 101
call ga_access(pars%g_erioo(isp),ilo,ihi,jlo,jhi,keri,ld)
call dtfsm('n','l','l','n','n',nri,jhi-jlo+1,1.0d0,w,
& dbl_mb(keri),nri)
call ga_release_update(pars%g_erioo(isp),1,nri,jlo,jhi)
101 continue
call ga_distribution(pars%g_erivv(isp),pars%me,ilo,ihi,jlo,jhi)
if (jlo.gt.jhi) cycle
call ga_access(pars%g_erivv(isp),ilo,ihi,jlo,jhi,keri,ld)
call dtfsm('n','l','l','n','n',nri,jhi-jlo+1,1.0d0,w,
& dbl_mb(keri),nri)
call ga_release_update(pars%g_erivv(isp),1,nri,jlo,jhi)
enddo
end

75
src/bsemol/bsemol.F Normal file
View file

@ -0,0 +1,75 @@
logical function bsemol(rtdb)
implicit none
#include "rtdb.fh"
integer rtdb
c
logical bse_doit
external bse doit
c
logical status
c
status = rtdb_cput(rtdb,'dft:theory', 1, 'dft')
if (status) status = bse_doit(rtdb)
bsemol = status
call grid_cleanup(.false.)
c
end function
c
c -----------------------------------------------------------------
c
logical function bse_doit(rtdb)
implicit none
#include "errquit.fh"
#include "global.fh"
#include "rtdb.fh"
#include "stdio.fh"
#include "util.fh"
#include "bse.fh"
integer, intent(in) :: rtdb
character(*), parameter :: pname = "bse_doit: "
type(bse_params_t) :: params
character*80 theory
integer me
c
c -----------------------------------------------------------------
c
c Initialize
c
bse_doit = .false.
c
c Determine the theory we are using
c
if (.not. rtdb_cget(rtdb, "dft:theory", 1, theory))
$ call errquit(pname//"theory not specified",0,RTDB_ERR)
c
c Initialize GW calculation
c
params%rtdb = rtdb
params%me = ga_nodeid()
c
call bse_init(params)
c
if (theory.eq."dft") then
if (params%davidson) then
call bse_davidson(params)
elseif (params%lanczos) then
call bse_lanczos(params)
else
call bse_analytic(params)
endif
c
elseif (theory.eq.'sodft') then
call errquit(pname//"spin-orbit not implemented",0,0)
else
call errquit(pname//"invalid theory: "//trim(theory),0,0)
endif
c
call bse_finalize(params)
bse_doit = .true.
c
c -----------------------------------------------------------------
c
end function

7
src/bsemol/scaleig.fh Normal file
View file

@ -0,0 +1,7 @@
c Copy of DFT scaleig
c
c $Id$
c
character*1 yev
integer nb
common /dft_scaleigp/ nb, yev

View file

@ -21,7 +21,7 @@ KNOWN_MODULE_SUBDIRS = NWints atomscf ddscf develop gradients moints nwdft \
property solvation nwpw fft rimp2_grad python argos analyz diana nwmd cafe space \
drdy vscf uccsdt qmmm rism qmd etrans tce geninterface bq mm cons\
perfm dntmc smd nbo dangchang leps ccca lucia dimqm rdmft fcidump\
gwmol xtb nwc_columbus
gwmol xtb nwc_columbus bsemol
# These are the libraries for the high-level modules. They should be
# specified in an order that will link correctly, but that shouldn't
@ -36,7 +36,7 @@ KNOWN_MODULE_LIBS = -ltest -lccsd -llucia -lmcscf -lselci -lmp2 -lmoints \
-lcafe -lspace -lanalyze -ldiana -lqhop \
-lpfft -ldplot -lnwpython -ldrdy -lvscf -lqmmm -lrism -lqmd -letrans\
-lpspw -ltce -lgeninterface -lbq -lmm -lcons -lperfm -lnbo -ldntmc -lccca\
-ldimqm -lrdmft -lfcidump -lgwmol -lnwxtb
-ldimqm -lrdmft -lfcidump -lgwmol -lnwxtb -lbsemol
# This include file handles configuration of the NW_MODULE_SUBDIRS and
# NW_MODULE_LIBS macros for what we actually want to build. It
@ -77,22 +77,22 @@ EXP_alldev = develop $(EXP_all)
EXP_pnnldev = develop $(EXP_pnnl)
EXP_qm = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \
optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans tce vscf dimqm \
gwmol
gwmol bsemol
EXP_smallqm = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \
optim stepper driver vib qmd mcscf property solvation dplot \
mp2_grad esp etrans vscf nwmd dimqm gwmol
mp2_grad esp etrans vscf nwmd dimqm gwmol bsemol
EXP_smallqmmm = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \
optim stepper driver vib qmd mcscf property solvation dplot \
mp2_grad esp etrans vscf nwmd qmmm prepar analyz dimqm gwmol
mp2_grad esp etrans vscf nwmd qmmm prepar analyz dimqm gwmol bsemol
EXP_notce = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \
optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans vscf dimqm \
gwmol
gwmol bsemol
EXP_tinyqmpw = nwdft property solvation driver ccsd mp2_grad hessian vib nwpw
EXP_ver2 = prepar gradients mp2_grad stepper vib mcscf \
optim driver analyz diana
#plane_wave ... broken? RJH
EXP_qmandpw = nwpw NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \
optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans vscf dimqm gwmol
optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans vscf dimqm gwmol bsemol
EXP_tps = prepar nwmd analyz diana
EXP_md = prepar nwmd analyz ddscf nwdft gradients esp
EXP_classical = prepar nwmd analyz
@ -187,6 +187,7 @@ MODULES := $(patsubst nwpw,nwpw fft,$(MODULES) )
MODULES := $(patsubst ddscf,ddscf atomscf nwdft NWints gradients,$(MODULES) )
MODULES := $(patsubst gwmol, gwmol nwdft ddscf gradients hessian atomscf NWints,$(MODULES) )
MODULES := $(patsubst bsemol, bsemol gwmol nwdft ddscf gradients hessian atomscf NWints,$(MODULES) )
MODULES := $(patsubst lucia, lucia mp2_grad ddscf nwdft property cphf gradients hessian atomscf NWints moints,$(MODULES) )
@ -274,6 +275,7 @@ MOD_LIBS := $(patsubst dimqm,-ldimqm,$(MOD_LIBS) )
MOD_LIBS := $(patsubst rdmft,-lrdmft,$(MOD_LIBS) )
MOD_LIBS := $(patsubst fcidump,-lfcidump,$(MOD_LIBS) )
MOD_LIBS := $(patsubst gwmol,-lgwmol,$(MOD_LIBS) )
MOD_LIBS := $(patsubst bsemol,-lbsemol,$(MOD_LIBS) )
MOD_LIBS := $(patsubst xtb,-lnwxtb,$(MOD_LIBS) )
MOD_LIBS := $(patsubst nwc_columbus,-lnwc_columbus,$(MOD_LIBS) )
MOD_LIBS := $(strip $(MOD_LIBS))
@ -339,6 +341,7 @@ STUBROUTINES := $(patsubst solvation, solvation,$(STUBROUTINES))
STUBROUTINES := $(patsubst dimqm,dimqm_addop dimqm_addop_uhf_damp dimqm_addop_uhf dimqm_tidy dimqm_input dimpar_input dimqm_main fock_dim dimqm_prep dim_grid_quadv0 dimqm_lclfld dimqm_polarizability dimqm_broadcast dimqm_Dens2Dipoles dimqm_rohf_calcDIM dimqm_rohf_wrap dimqm_used dimqm_setgo dimqm_getgo dimqm_getenergy dimqm_setlfirst dimqm_getlrsp dimqm_lrsp dimqm_getlclfld dimqm_getgdipeli dimqm_energy dimqm_init,$(STUBROUTINES) )
#STUBROUTINES := $(patsubst cons, cons_end cons_input cons_add_egrad cons_add_energy cons_init cons_print_summary cons_active,$(STUBROUTINES) )
STUBROUTINES := $(patsubst gwmol, gw_input gwmol, $(STUBROUTINES) )
STUBROUTINES := $(patsubst bsemol, bse_input bsemol, $(STUBROUTINES) )
STUBROUTINES := $(patsubst xtb, xtb_input xtb_energy xtb_energy_gradient, $(STUBROUTINES) )
STUBROUTINES := $(patsubst nwc_columbus, sifcfg sifwh sodft_aoints, $(STUBROUTINES) )
STUBROUTINES := $(strip $(STUBROUTINES))

View file

@ -308,16 +308,7 @@ c
write(LuOut,*) ' B matrix '
call output(x, 1, ncart, 1, nzvar, ncart, nzvar, 1)
endif
call dgemm('t', 'n', nzvar, nzvar, ncart, 1d0,
$ x, ncart, x, ncart, 0d0, y, nzvar) ! BT.B -> y ... should symmetrize
do i = 1, nzvar
do j = 1, nzvar
ij = i + (j-1)*nzvar
ji = j + (i-1)*nzvar
y(ij) = 0.5d0*(y(ij) + y(ji))
y(ji) = y(ij)
end do
end do
call dsyrk('u', 't', nzvar, ncart, 1d0, x, ncart, 0d0, y, nzvar)
call dsyev('v', 'u', nzvar, y, nzvar, e, z, ndim**2, info)
if (info .ne. 0) call errquit('geom_binvr: dsyev failed', info,
& MEM_ERR)

View file

@ -12,6 +12,7 @@ c
type(gw_params_t) :: pars
c
character(len=13), parameter :: pname = 'gw_analytic: '
character*255 gwevals
c
integer,dimension(2) :: nocc,nvir,npoles,nstates,mynpoles
@ -45,6 +46,10 @@ c
integer maxpoles,maxnpoles,isp,lSize
logical ok
integer,parameter :: unitno = 66
logical,external :: movecs_write
c
c -----------------------------------------------------------------
c
@ -89,6 +94,15 @@ c
call gw_vxc(dbl_mb(k_sigma),nmo,pars)
if (pars%ncap) then
do isp=1,ipol
call daxpy(pars%nocc(isp),1d0,pars%vxddp(isp),0,
$ dbl_mb(k_sigma+(isp-1)*pars%nmo),1)
call daxpy(pars%nvir(isp),1d0,pars%vxddm(isp),0,
$ dbl_mb(k_sigma+(isp-1)*pars%nmo+pars%nocc(isp)),1)
enddo
endif
pars%timer%vxc = pars%timer%vxc + util_wallsec() - timer
c
c Obtain Sigma_x
@ -135,6 +149,14 @@ c
endif
write(luout,*)
endif
call gw_shift(dbl_mb(k_eold),dbl_mb(pars%k_evals),-efermi,
$ pars%nmo,nbf_ao,pars%ipol)
if(.not.movecs_write(pars%rtdb, ao_bas_han, movecs_out, 'dft',
$ 'W vals',nbf_ao,ipol,(/pars%nmo,pars%nmo/),
$ dbl_mb(pars%k_occ), nbf_ao,
$ dbl_mb(pars%k_evals),nbf_ao,pars%g_movecs))
$ call errquit('gw_analytic: could not save movecs',0,DISK_ERR)
c
c Obtain eigenvalue difference
c
@ -185,25 +207,26 @@ c
endif
c
if (pars%evgw0.or.pars%evgw) then
call gw_scissor(dbl_mb(k_Enew),dbl_mb(k_Eold),pars%nmo,
call gw_scissor(dbl_mb(k_Enew),dbl_mb(k_Eold),pars%nmo,
$ pars%nocc,pars%nvir,pars%noqp,pars%nvqp,
$ pars%ipol,pars%me)
call ga_brdcst(1038,dbl_mb(k_enew),
$ ma_sizeof(mt_dbl,nmo,mt_byte),0)
call ga_brdcst(1038,dbl_mb(k_sigma),
$ ma_sizeof(mt_dbl,nmo,mt_byte),0)
call ga_brdcst(1038,dbl_mb(k_enew),
$ ma_sizeof(mt_dbl,pars%ipol*nmo,mt_byte),0)
call ga_brdcst(1038,dbl_mb(k_sigma),
$ ma_sizeof(mt_dbl,pars%ipol*nmo,mt_byte),0)
if (pars%evgw0.or.pars%evgw) then
!
! Update Fermi energy
!
do i=1,ipol
lSize = k_enew + (i-1)*nbf_ao
lSize = k_enew + (i-1)*nmo
efermi_new(i) = gw_efermi(dbl_mb(lSize),pars%nocc(i))
enddo
call gw_shift(dbl_mb(k_enew),dbl_mb(k_eold),efermi_new,
$ nbf_ao,pars%nmo,pars%ipol)
$ pars%nmo,pars%nmo,pars%ipol)
efermi(:) = efermi(:) + efermi_new(:)
call dcopy(pars%ipol*nmo,dbl_mb(k_eold),1,dbl_mb(k_enew),1)
@ -214,6 +237,23 @@ c
endif
endif
!
! Save eigenvalues
!
if (main) then
call gw_shift(dbl_mb(k_enew),dbl_mb(k_eold),-efermi,nmo,nmo,ipol)
call util_file_name('gwevals',.false.,.false.,gwevals)
call util_file_name_resolve(gwevals,.false.)
open(unit=unitno,status='unknown',form='unformatted',file=gwevals)
do i=1,ipol
call swrite(unitno, dbl_mb(k_eold+(i-1)*nmo),nmo)
enddo
close(unit=unitno)
endif
9050 format(10x,I3,6x,F8.3,6x,F6.3,6x,F8.3)
c
if (.not.ma_chop_stack(l_sigma))

View file

@ -33,7 +33,7 @@
integer nclusters, clusters(pars%nmo), myll, myul
integer llqp, ulqp, icluster
double precision upper, lower, resupper, reslower, const
double precision values(20), errors(20), radius, maxradius
double precision values(100), errors(100), radius, maxradius
c
!
! Initialize variables
@ -60,7 +60,7 @@ c !
if (eviter.gt.4) then
maxradius = 0.010d0
else
maxradius = 0.1d0
maxradius = 0.030d0
endif
c !
@ -75,9 +75,6 @@ c !
totalqp = pars%noqp(isp) + pars%nvqp(isp)
if (totalqp.lt.1) cycle
totalqp = pars%noqp(isp) + pars%nvqp(isp)
if (totalqp.lt.1) cycle
esterror(:) = 0d0
if ((isp.eq.1).and.main) then
@ -139,16 +136,18 @@ c !
if (eviter.lt.10) then
if (myul.lt.ulqp) then
Eout = Enew(llmo+iqp,isp) + 0.002d0
c Eout = Enew(llmo+iqp,isp) + 0.002d0
Eout = Eout + 0.001d0
elseif (myll.gt.llqp) then
Eout = Enew(llmo+iqp-2,isp) - 0.002d0
c Eout = Enew(llmo+iqp-2,isp) - 0.002d0
Eout = Eout - 0.001d0
endif
endif
bracket = .false.
const = Eold(llmo+iqp-1,isp) + Sigma_X - sigmaold
lower = 1d100
upper = -1d100
lower = -1d100
upper = 1d100
c Compute Screened-Coulomb matrix elements
timer = util_wallsec()
@ -164,7 +163,7 @@ c Compute Galiitski-Migdal correlation energy
endif
c
do iter=1,pars%qpiter+5
do iter=1,pars%qpiter+15
timer2 = util_wallsec()
Ein = Eout

View file

@ -48,6 +48,10 @@ c
logical converged,ok
double precision, external :: gw_efermi
logical,external :: movecs_write
integer,parameter :: unitno = 66
character*255 gwevals
c
c -----------------------------------------------------------------
c
@ -83,6 +87,15 @@ c
call gw_vxc(dbl_mb(kSigma),nmo,pars)
if (pars%ncap) then
do isp=1,ipol
call daxpy(pars%nocc(isp),1d0,pars%vxddp(isp),0,
$ dbl_mb(ksigma+(isp-1)*pars%nmo),1)
call daxpy(pars%nvir(isp),1d0,pars%vxddm(isp),0,
$ dbl_mb(ksigma+(isp-1)*pars%nmo+pars%nocc(isp)),1)
enddo
endif
pars%timer%vxc = pars%timer%vxc + util_wallsec() - timer
c
c Obtain Sigma_x
@ -127,7 +140,6 @@ c
lSize = lSize*max(pars%eviter,1)
if (.not.ma_push_get(mt_dbl,lSize,'values',lValues,kValues))
$ call errquit(pname//'failed to allocate values',0,MA_ERR)
lSize = lSize*max(pars%eviter,1)
if (.not.ma_push_get(mt_dbl,lSize,'errors',lErrors,kErrors))
$ call errquit(pname//'failed to allocate errors',0,MA_ERR)
endif
@ -169,19 +181,14 @@ c
9020 format(/,2x,"*** G",I1,"W",I1)
9021 format(/,2x,"*** G",I2,"W",I2)
if (eviter.eq.0 .and. pars%ncap .and. pars%evgw) then
do i=1,ipol
do j=0,pars%nmo-1
if (j.lt.pars%nocc(i)) then
dbl_mb(kEmod+(i-1)*nmo+j) =
$ dbl_mb(kEmod+(i-1)*nmo+j) + pars%vxddp(i)
else
dbl_mb(kEmod+(i-1)*nmo+j) =
$ dbl_mb(kEmod+(i-1)*nmo+j) + pars%vxddm(i)
endif
enddo
enddo
endif
call gw_shift(dbl_mb(kemod),dbl_mb(pars%k_evals),-efermi,
$ pars%nmo,nbf_ao,pars%ipol)
if(.not.movecs_write(pars%rtdb, ao_bas_han, movecs_out, 'dft',
$ 'W vals',nbf_ao,ipol,(/pars%nmo,pars%nmo/),
$ dbl_mb(pars%k_occ), nbf_ao,
$ dbl_mb(pars%k_evals),nbf_ao,pars%g_movecs))
$ call errquit('gw_analytic: could not save movecs',0,DISK_ERR)
c
c Obtain eigenvalue difference
c
@ -194,23 +201,25 @@ c
c
c Experimental bit, might be useful for BSE calculations
c
if (eviter.eq.0 .and. pars%ncap .and. (.not.pars%evgw)) then
do i=1,ipol
do j=0,pars%nmo-1
if (j.lt.pars%nocc(i)) then
dbl_mb(kEmod+(i-1)*nmo+j) =
$ dbl_mb(kEmod+(i-1)*nmo+j) + pars%vxddp(i)
else
dbl_mb(kEmod+(i-1)*nmo+j) =
$ dbl_mb(kEmod+(i-1)*nmo+j) + pars%vxddm(i)
endif
enddo
enddo
endif
c if (eviter.eq.0 .and. pars%ncap .and. (.not.pars%evgw)) then
c do i=1,ipol
c do j=0,pars%nmo-1
c if (j.lt.pars%nocc(i)) then
c dbl_mb(kEmod+(i-1)*nmo+j) =
c $ dbl_mb(kEmod+(i-1)*nmo+j) + pars%vxddp(i)
c else
c dbl_mb(kEmod+(i-1)*nmo+j) =
c $ dbl_mb(kEmod+(i-1)*nmo+j) + pars%vxddm(i)
c endif
c enddo
c enddo
c endif
c
c if (eviter.eq.0 .and. pars%ncap) then
c call dcopy(nmo*ipol,dbl_mb(kEmod),1,dbl_mb(kEref),1)
c endif
c
call dcopy(nmo*ipol, dbl_mb(kEmod), 1, dbl_mb(kEref), 1)
c
c Obtain Screened Coulomb matrices on imaginary axis
c
@ -277,11 +286,11 @@ c
! Update Fermi energy
!
do i=1,ipol
lSize = kEnew + (i-1)*nbf_ao
lSize = kEnew + (i-1)*nmo
efermi_new(i) = gw_efermi(dbl_mb(lSize),pars%nocc(i))
enddo
call gw_shift(dbl_mb(kEnew),dbl_mb(kemod),efermi_new,
$ nbf_ao,pars%nmo,pars%ipol)
$ pars%nmo,pars%nmo,pars%ipol)
efermi(:) = efermi(:) + efermi_new(:)
call dcopy(pars%ipol*nmo,dbl_mb(kemod),1,dbl_mb(kEnew),1)
@ -293,12 +302,26 @@ c
endif
endif
if (main) then
call gw_shift(dbl_mb(kenew),dbl_mb(kemod),-efermi,nmo,nmo,ipol)
call util_file_name('gwevals',.false.,.false.,gwevals)
call util_file_name_resolve(gwevals,.false.)
open(unit=unitno,status='unknown',form='unformatted',file=gwevals)
do i=1,ipol
call swrite(unitno, dbl_mb(kemod+(i-1)*nmo),nmo)
enddo
close(unit=unitno)
endif
if (.not.ga_destroy(pars%g_sols))
$ call errquit('gw_cdgw_newton: cannot destroy GA',0,GA_ERR)
if (.not.(ma_free_heap(liW).and.ma_free_heap(liWidx)))
$ call errquit(pname//'cannot destroy iW arrays',0,MA_ERR)
c
end subroutine

View file

@ -60,7 +60,7 @@ C
logical main, ok, parallel_diag, skip
logical, external :: is_near
double precision, parameter :: shift = 0.0001d0
double precision, parameter :: shift = 0.00005d0
double complex :: w, R1, dR1
double complex :: ieta

View file

@ -40,7 +40,7 @@
integer nclusters, clusters(pars%nmo),myll,myul
integer llqp,ulqp,icluster
double precision upper,lower,resupper,reslower,const
double precision values(25),errors(25),radius
double precision values(100),errors(100),radius
main = pars%me .eq. 0
@ -110,9 +110,11 @@
if (eviter.lt.3) then
if (myul.lt.ulqp) then
Eout = Enew(llmo+iqp,isp) + 0.001d0
!Eout = Enew(llmo+iqp,isp) + 0.001d0
Eout = Eout + 0.001d0
elseif (myll.gt.llqp) then
Eout = Enew(llmo+iqp-2,isp) - 0.001d0
!Eout = Enew(llmo+iqp-2,isp) - 0.001d0
Eout = Eout - 0.001d0
endif
endif
@ -207,6 +209,7 @@ c
! Save last energy
Enew(llmo+iqp-1,isp) = Eout
Sigma(llmo+iqp-1,isp) = SigmaC + Sigma_X
! check convergence info
if (converged) then

View file

@ -27,21 +27,17 @@ c
& rtdb_put(rtdb,'gw:graph',mt_log,1,.false.).and.
& rtdb_put(rtdb,'gw:linear',mt_log,1,.false.).and.
& rtdb_put(rtdb,'gw:eviter',mt_int,1,4).and.
& rtdb_put(rtdb,'gw:qpiter',mt_int,1,10).and.
& rtdb_put(rtdb,'gw:qpiter',mt_int,1,20).and.
& rtdb_put(rtdb,'gw:noqp_alpha',mt_int,1,1).and.
& rtdb_put(rtdb,'gw:nvqp_alpha',mt_int,1,0).and.
& rtdb_put(rtdb,'gw:noqp_beta',mt_int,1,1).and.
& rtdb_put(rtdb,'gw:nvqp_beta',mt_int,1,0).and.
& rtdb_put(rtdb,'gw:ngl',mt_int,1,200).and.
& rtdb_put(rtdb,'gw:ngl',mt_int,1,160).and.
& rtdb_put(rtdb,'gw:nomega',mt_int,1,100).and.
& rtdb_put(rtdb,'gw:first',mt_int,1,1).and.
& rtdb_put(rtdb,'gw:threshold',mt_dbl,1,0.005d0/ha2ev).and.
& rtdb_put(rtdb,'gw:ncap',mt_log,1,.false.).and.
#ifdef GWCMPLX
& rtdb_put(rtdb,'gw:eta',mt_dbl,1,0.0001d0).and.
#else
& rtdb_put(rtdb,'gw:eta',mt_dbl,1,0.00001d0).and.
#endif
& rtdb_put(rtdb,'gw:eta',mt_dbl,1,0.01d0).and.
& rtdb_put(rtdb,'gw:rpa',mt_log,1,.false.).and.
& rtdb_put(rtdb,'gw:diag',mt_log,1,.false.).and.
& rtdb_put(rtdb,'gw:core',mt_log,1,.false.).and.

View file

@ -30,8 +30,12 @@ C
!
!take full newton step for large enough z
!
if (z.gt.0.4d0 .and. z.le.1d0 ) then
xout = xin + step
if (z.gt.0.6d0 .and. z.le.1d0 ) then
if (iter.eq.1) then
xout = xin + 0.2d0*step
else
xout = xin + 0.7d0*step
endif
!
! when we have bracketed the solution, the
@ -51,7 +55,11 @@ C
! Take an scaled Newton step for small Z
!
elseif (z.gt.0.1d0) then
xout = xin + 0.60d0*step
if (iter.eq.1) then
xout = xin + 0.10d0*step
else
xout = xin + 0.30d0*step
endif
!
! If Z is greater than 1.0 or very small we have to suspect
@ -59,7 +67,8 @@ C
! steepest decent step
!
else
xout = xin + radius*fx
c xout = xin + radius*fx
xout = xin + sign(0.002d0,fx)
endif

View file

@ -31,6 +31,8 @@ c
c
logical movecs_read, movecs_read_header, atom_tag_check
external movecs_read, movecs_read_header, atom_tag_check
logical,external :: movecs_write
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
@ -46,8 +48,6 @@ c
c
c -----------------------------------------------------------------
c
zeta = 0.304121d0
aq2 = 2d0/9d0 * (1d0-zeta)**2
pars%timer%init = util_wallsec()
pars%timer%total = pars%timer%init
@ -71,6 +71,9 @@ c
pars%me = ga_nodeid()
pars%nprocs = ga_nnodes()
!oskel = .false.
!oadapt = .false.
#ifdef USE_OPENMP
pars%iMaxThreads = omp_get_max_threads()
#else
@ -135,9 +138,9 @@ c
$ pars%nvqp(1)))
$ call errquit(pname//'failed to read nvqp',0,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb,'gw:noqp_beta',mt_int,1,pars%noqp(2)))
$ call errquit(pname//'failed to read noqp',0,RTDB_ERR)
$ call errquit(pname//'failed to read noqp beta',0,RTDB_ERR)
if (.not.rtdb_get(pars%rtdb,'gw:nvqp_beta',mt_int,1,pars%nvqp(2)))
$ call errquit(pname//'failed to read nvqp',0,RTDB_ERR)
$ call errquit(pname//'failed to read nvqp beta',0,RTDB_ERR)
c
if (.not.rtdb_get(pars%rtdb,'gw:analyticw',mt_log,1,pars%anaw))
$ call errquit(pname//'failed to read analyticw',0,RTDB_ERR)
@ -236,7 +239,6 @@ c
do j=1,pars%nmo
if (dbl_mb(pars%k_occ+(i-1)*nbf_ao+j-1).gt.0.d0) then
nocc = nocc + 1
cycle
else
nvir = nvir + 1
endif
@ -244,7 +246,7 @@ c
pars%nocc(i) = nocc
pars%nvir(i) = nvir
if (nocc+nvir.ne.pars%nmo)
$ call errquit(pname//'something went wrong',0,0)
$ call errquit(pname//'Nocc + Nvir /= NMO ',0,0)
c
if ((pars%noqp(i).gt.nocc).or.(pars%noqp(i).eq.-1))
@ -253,12 +255,11 @@ c
if (pars%noqp(i).lt.-1)
$ call errquit(pname//'invalid number noqp specified',0,0)
c
if ((pars%nvqp(i).gt.nvir).or.(pars%nvqp(i).eq.-1))
if ((pars%nvqp(i).gt.nvir).or.(pars%nvqp(i).le.-1))
$ pars%nvqp(i) = nvir
if (pars%nvqp(i).lt.-1)
$ call errquit(pname//'invalid number nvqp specified',0,0)
! If core option was passed, then start counting from the bottom
! else, start counting from the top (Fermi level)
if (pars%core) then
pars%llmo(i) = pars%first
else
@ -284,16 +285,23 @@ c
c
1001 format(2X,A12,F8.3,A3)
if (pars%ncap) then
zeta = 0.304121d0
if ((xfac(87)).gt.1d-10) zeta=0.5d0
aq2 = 2d0/9d0 * (1d0-zeta)**2
do i=1,ipol
ehomo = dbl_mb(pars%k_evals+(i-1)*nbf_ao+pars%nocc(i)-1)
vxddp = -0.5d0*aq2*(1d0 + sqrt(1d0-4d0*ehomo/aq2))
vxddm = -0.5d0*aq2*(1d0 - sqrt(1d0-4d0*ehomo/aq2))
vxddp = -0.5d0*aq2*(1d0 + sqrt(1d0-4d0*ehomo/aq2))!*0.15d0
vxddm = -0.5d0*aq2*(1d0 - sqrt(1d0-4d0*ehomo/aq2))!*0.15d0
pars%vxddp(i) = vxddp
pars%vxddm(i) = vxddm
if (pars%me.eq.0) then
write(luout,1001) ' v_x^{DD+}: ',vxddp*ha2ev, ' eV'
write(luout,1001) ' v_x^(DD-): ',vxddm*ha2ev, ' eV'
endif
call daxpy(pars%nocc(i),1d0,vxddp,0,
$ dbl_mb(pars%k_evals+(i-1)*nbf_ao),1)
call daxpy(pars%nvir(i),1d0,vxddm,0,
$ dbl_mb(pars%k_evals+(i-1)*nbf_ao+pars%nocc(i)),1)
enddo
endif
c

View file

@ -484,6 +484,7 @@ c if (.not.ga_create(mt_dbl,nri,npairs,erilabel,nri,0,eri))
enddo
ok = ma_chop_stack(l_tmp)
ok = ma_chop_stack(lmap)
call ga_distribution(eri,me,ld,ld,polelo,polehi)

View file

@ -93,6 +93,7 @@ c
$ pars%g_movecs(isp),pars%g_movecs(isp),0d0,
$ g_dens(isp))
if (odftps) call pstat_off(ps_dgemm)
call ga_symmetrize(g_dens(isp))
enddo
c
nscr = 2*nbf_ao_mxnbf_ce**2
@ -131,7 +132,7 @@ c
if (.not.ma_push_get(mt_dbl,2*nbf_ao,'vxc temp',l_tmp,k_tmp))
$ call errquit(pname//'failed to allocate temp',0,MA_ERR)
do isp=1,ipol
if (oskel)
if (oskel)
$ call sym_symmetrize(geom, ao_bas_han, .false., g_vxc(isp))
call ga_dgemm('n','n',nbf_ao,nmo,nbf_ao,1.0d0,g_vxc(isp),
$ pars%g_movecs(isp),0.0d0,g_dens(isp))

View file

@ -332,6 +332,8 @@ ckbn mrcc-r-1 -3
#endif
else if (inp_compare(.false.,test,'gw')) then
call gw_input(rtdb)
else if (inp_compare(.false.,test,'bse')) then
call bse_input(rtdb)
else if (inp_compare(.false.,test,'xtb')) then
call xtb_input(rtdb)
c

View file

@ -703,6 +703,7 @@ c
xname(83)=' MN15-L Exchange Functional'
xname(84)=' MN15 Exchange Functional'
xname(85)=' M06-SX Exchange Functional'
xname(87)=' NCAPR GGA Exchange Functional'
c
c Exchange-Correlation combination functional name defaults.
c
@ -785,6 +786,7 @@ c
xcname(84)=' MN15 Method XC Functional'
xcname(85)=' M06-SX Method XC Functional'
xcname(86)=' r^2SCAN0 Method XC Functional'
xcname(87)=' NCAPR Method XC Functional'
c
c place character XC parameters in rtdb
c

View file

@ -153,7 +153,7 @@ c
c r2scan0 ==> 19
c
integer num_dirs, ind, mlen, iline, n, itemp
parameter (num_dirs=168)
parameter (num_dirs=169)
c
character*15 dirs(num_dirs)
character*255 test
@ -212,7 +212,7 @@ c
& 'xscan', 'cscan',
& 'hcth147@tz2p',
W 'wb97','wb97x','wb97x-d3',
& 'xr2scan', 'cr2scan','r2scan0' /
& 'xr2scan', 'cr2scan','r2scan0', 'ncapr' /
c becke97-d, ssb-d
logical disp, dumdirect
integer ivdw
@ -298,7 +298,8 @@ c
& 9710,9720,9650,9660,9615,9616,
& 1947,
& 9410,9420,9430,
& 9661,9662,9670,
& 9661,9662,9670,
& 9730,
& 1999 ) ind
@ -1817,6 +1818,17 @@ c
nlxfac(72) = .true.
xfac(72) = 1.0d0
goto 10
9730 xccomb(87) = .true.
cfac(3) = 1d0
lcfac(3) = .true.
cfac(4) = 1d0
nlcfac(4) = .true.
lxfac(87) = .true.
nlxfac(87) = .true.
xfac(87) = 1.0d0
goto 10
c
c r2SCAN metaGGA
c

View file

@ -1646,6 +1646,7 @@ c Clean scratch space
if (nodezero) call util_flush(LuOut)
c
endif ! gmh
c
enddo ! outer loop over roots
c

View file

@ -1855,15 +1855,32 @@ c
if (.not. do_2nd .and..not. do_3rd) then
call xc_xncap(tol_rho, xfac(72), rho,
& delrho, Amat, Cmat, nq, ipol,
& Ex, qwght, ldew, func)
& Ex, qwght, ldew, func, ' ')
else if (.not. do_3rd) then
call xc_xncap_d2(tol_rho, xfac(72), rho,
& delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol,
& Ex, qwght, ldew, func)
& Ex, qwght, ldew, func, ' ')
else
call xc_xncap_d3(tol_rho, xfac(72), rho,
& delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3,
& nq, ipol, Ex, qwght, ldew, func)
& nq, ipol, Ex, qwght, ldew, func, ' ')
endif
endif
c
if (abs(xfac(87)).gt.eps)then
if (.not. do_2nd .and..not. do_3rd) then
call xc_xncap(tol_rho, xfac(87), rho,
& delrho, Amat, Cmat, nq, ipol,
& Ex, qwght, ldew, func, 'r')
else if (.not. do_3rd) then
call xc_xncap_d2(tol_rho, xfac(87), rho,
& delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol,
& Ex, qwght, ldew, func, 'r')
else
call xc_xncap_d3(tol_rho, xfac(87), rho,
& delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3,
& nq, ipol, Ex, qwght, ldew, func, 'r')
endif
endif
c

View file

@ -32,19 +32,23 @@ c
c
#if !defined SECOND_DERIV && !defined THIRD_DERIV
Subroutine xc_xncap(tol_rho, fac, rho, delrho,
& Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
& Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,
& version)
#elif defined(SECOND_DERIV) && !defined THIRD_DERIV
Subroutine xc_xncap_d2(tol_rho, fac, rho, delrho,
& Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
& qwght,ldew,func)
& qwght,ldew,func,
& version)
#else
Subroutine xc_xncap_d3(tol_rho, fac, rho, delrho,
& Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3,
& nq, ipol, Ex, qwght,ldew,func)
& nq, ipol, Ex, qwght,ldew,func,
& version)
#endif
c
implicit none
c
c
character*1 version
double precision fac, Ex
integer nq, ipol
logical ldew
@ -77,8 +81,7 @@ c
double precision ckf, Ax
double precision F43, F13, F23, F49
c
parameter(mu=0.2195149727645171d0, beta=0.01808569669d0)
parameter(zeta=0.30412141859531383d0)
parameter(mu=0.2195149727645171d0)
parameter (F43=4.d0/3.d0, F13=1.d0/3.d0, F23=2.0d0/3.0d0)
parameter (F49=4d0/9d0)
@ -100,11 +103,22 @@ c
#endif
double precision ops, omz, logops
double precision s2overgamma,s2overgamma_rho43
double precision dFxd1s,dFxd2s
arcsinh(s)=dlog(s+dsqrt(1d0+s*s))
darcsinh(s)=1d0/dsqrt(1d0+s*s)
d2arcsinh(s) = -s/dsqrt(1d0+s*s)**3
d3arcsinh(s) = (2d0*s*s - 1d0)/dsqrt(1d0+s*s)**5
c
if (version.eq.'r') then
zeta = 0.5d0
beta = 0.017983d0
else
zeta = 0.30412141859531383d0
beta = 0.01808569669d0
endif
pi = acos(-1.d0)
ckf = (3d0*pi*pi)**F13
Ax = -3d0/(4d0*pi)*ckf
@ -130,11 +144,21 @@ CDEC$ NOSWP
s = dsqrt(gamma)/(2d0*ckf*rho43)
s2 = s*s
s2overgamma = 1d0/(4d0*ckf*ckf*rho43*rho43)
s2overgamma_rho43 = 1d0/(4d0*ckf*ckf*rho43)
d1s = 0.5d0*s/gamma
if (s.lt.1d-8) then
Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s +
& mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2
dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 +
& 4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s
dFxd1s = mu*s2overgamma*
& (1d0 + 1.5d0*alpha*zeta*s +
& 2d0*(alpha-beta-alpha*zeta-0.5d0)*s2)
dFxd2s = -0.5d0*mu*
& (1d0 + 1.5d0*alpha*zeta*s +
& 2d0*(alpha-beta-alpha*zeta-0.5d0)*s2)
d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s +
& 12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2
d3Fx = 6d0*alpha*zeta*mu +
@ -147,7 +171,7 @@ CDEC$ NOSWP
asymp = 1d0 + alpha*(omz*s*logops + zeta*s)
dasymp = alpha*( (s + zeta)/ops + omz*logops )
if (s.gt.30d0) then
if (s.gt.20d0) then
tansin = arcsinh(s)
dtansin = darcsinh(s)
else
@ -159,18 +183,19 @@ CDEC$ NOSWP
ddenom = beta*dtansin
factor = tansin/denom
dfactor = (dtansin - factor*ddenom)/denom
dfactor = dtansin*(1d0 - beta*factor)/denom
Fx = 1d0 + mu*factor*asymp
dFx = mu*(dfactor*asymp + factor*dasymp)
dFxd1s = dFx*d1s
dFxd2s = -0.5d0*dFxd1s/s2overgamma
#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
if (s.gt.30d0) then
if (s.gt.24d0) then
d2tansin = d2arcsinh(s)
else
d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 +
& dtanh(s)*d2arcsinh(s) -
& 2d0*tansin/dcosh(s)**2
d2tansin = (darcsinh(s)-tansin)*2d0/dcosh(s)**2 +
& dtanh(s)*d2arcsinh(s)
endif
d2asymp = alpha*omz*(2d0 + s)/ops**2
d2denom = beta*d2tansin
@ -202,9 +227,8 @@ CDEC$ NOSWP
Ex = Ex + Ax*rho43*Fx*qwght(n)*fac
if (ldew) func(n) = func(n) + Ax*rho43*Fx*fac
d1s = 0.5d0*s/gamma
Amat(n,1) = Amat(n,1) + F43*Ax*rho13*(Fx - s*dFx)*fac
Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*Ax*rho43*dFx*d1s*fac
Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*Ax*rho43*dFxd1s*fac
#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
@ -214,10 +238,10 @@ CDEC$ NOSWP
& (Fx - s*dFx + 4d0*s2*d2Fx)*fac
Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) -
& (F43*Ax*rho13*d2Fx*d1s*s)*4d0*fac
& (0.5d0*F43*Ax*rho13*d2Fx*s2overgamma)*4d0*fac
Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + Ax*rho43*
& (d2Fx*d1s**2 + dFx*d2s)*8d0*fac
Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + Ax*
& s2overgamma_rho43/gamma*(0.25d0*d2Fx + dFxd2s)*8d0*fac
#endif
#ifdef THIRD_DERIV
rhom53 = rhom23*rrho

View file

@ -9,7 +9,7 @@
task_save_state.o task.o task_input.o task_shell_input.o task_num_grad.o \
task_vscf.o task_property.o task_dntmc.o task_bsse.o \
task_jefftce.o task_ncc.o task_rfcf.o \
task_gwmol.o \
task_gwmol.o task_bsemol.o \
task_aoints.o task_nadtrc.o
USES_BLAS = task_num_grad.F task_bsse.F

View file

@ -163,6 +163,10 @@ c
logical task_nadtrc
external task_nadtrc
logical task_bsemol
external task_bsemol
c logical embed_optimize
c external embed_optimize
c
@ -468,6 +472,8 @@ c status=embed_prepare(rtdb)
status=task_dntmc(rtdb)
else if (inp_compare(.false., 'gw', operation)) then
status = task_gwmol(rtdb)
else if (inp_compare(.false., 'bse', operation)) then
status = task_bsemol(rtdb)
c
c **** check for theories that have specific operations ****
c

138
src/task/task_bsemol.F Normal file
View file

@ -0,0 +1,138 @@
c
c Task handler for Molecular BSE calcualtions
c Daniel Mejia-Rodriguez, 2022
c
logical function task_bsemol(rtdb)
c
implicit none
c
#include "errquit.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "inp.fh"
#include "stdio.fh"
#include "util.fh"
#include "global.fh"
c
integer rtdb
c
character*32 theory
character*255 title
integer iproc, ipol
logical oprint, rodft, dftmp2, status, converged, ignore
double precision cputime, walltime
c
logical bsemol, task_energy
external bsemol, task_energy
c
c Initialization
c
task_bsemol = .false.
c
c Time at the start
c
cputime = util_cpusec()
walltime = util_wallsec()
c
c Invalidate status
c
if (.not. rtdb_put( rtdb, 'task:status', mt_log, 1, .false.))
& call errquit('task_bse: failed to invalidate status',0,RTDB_ERR)
c
c Print BSE header
c
oprint = util_print('information', print_low)
if (oprint) then
if (.not. rtdb_cget(rtdb, 'title', 1, title)) title = ' '
if (ga_nodeid().eq.0) then
write(LuOut,*)
call util_print_centered(
& LuOut, 'NWChem BSE Module', 40, .true.)
write(LuOut,*)
write(LuOut,*)
if (title.ne.' ') then
call util_print_centered(LuOut, title, 40, .false.)
write(LuOut,*)
write(LuOut,*)
endif
endif
endif
c
c Check if theory has been specified
c
if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory))
& call errquit('task_bsemol: theory not specified',0, RTDB_ERR)
c
c Restrict to DFT for now (SCF could also be used)
c
if (.not.inp_compare(.false.,'dft',theory))
& call errquit('task_bsemol: BSE is only implemented for DFT', 0,
& INPUT_ERR)
c
if (rtdb_get(rtdb, 'dft:rodft', mt_log, 1, rodft))
& call errquit('task_bsemol: BSE is not implemented for RODFT', 0,
& INPUT_ERR)
c
if (rtdb_get(rtdb, 'dft:dftmp2', mt_log, 1, dftmp2))
& call errquit('task_bsemol: BSE is not implemented for double'//
& ' hybrid functionals',0,INPUT_ERR)
c
if (.not.rtdb_get(rtdb, 'dft:ipol', mt_int, 1, ipol))
& call errquit('task_bsemol: failed to read ipol', 0, RTDB_ERR)
if (.not. rtdb_get(rtdb, 'task:ignore', mt_log, 1, ignore))
& ignore = .true.
c
c Make sure a DFT solution is available
c
status = rtdb_get(rtdb, 'dft:converged', mt_log, 1, converged)
if (.not.status) then
status = task_energy(rtdb)
if (.not.rtdb_get(rtdb, 'dft:converged', mt_log, 1, converged))
& call errquit('task_bsemol: failed to read convergence info',0,
& RTDB_ERR)
endif
status = converged
if (ga_nodeid().eq.0) then
write(luout,*) " "
write(luout,'(29X,A)') "***********************"
write(luout,'(29X,A)') "*** BSE ***"
write(luout,'(29X,A)') "***********************"
write(luout,*) " "
if (.not.converged) then
write(LuOut,*)
write(LuOut,'(27X,A)') '!!! WARNING: !!!'
write(LuOut,'(27X,A)') '!!! DFT DID NOT CONVERGE !!!'
write(LuOut,*)
endif
call util_flush(luout)
endif
if ((.not.status) .and. (.not.ignore))
& call errquit('task_bsemol: Energy failure',0,CALC_ERR)
c
c Actually compute GW
c
status = bsemol(rtdb)
c
cputime = util_cpusec() - cputime
walltime = util_wallsec() - walltime
c
if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .true.))
& call errquit('task_bse: failed to set status',0,
& RTDB_ERR)
if (.not. rtdb_put(rtdb, 'task:cputime', mt_dbl, 1, cputime))
& call errquit('task_bse: failed storing cputime',0,
& RTDB_ERR)
if (.not. rtdb_put(rtdb, 'task:walltime', mt_dbl, 1, walltime))
& call errquit('task_bse: failed storing walltime',0,
& RTDB_ERR)
c
task_bsemol = status
c
end function

View file

@ -37,7 +37,7 @@ c
c task [qmmm] <string theory> [<string operation = energy>] [numerical || analytic]
c [ignore]
c
logical numerical, ignore, oqmmm, analytic
logical numerical, ignore, oqmmm, analytic, gwbse
logical orism
logical omm
logical opes
@ -61,6 +61,7 @@ c
opes = .false.
ofragment = .false.
osmd = .false.
gwbse = .false.
c oembed = .false.
c
call rtdb_clean(rtdb, 'task:') ! Delete ALL previous task parameters
@ -114,6 +115,10 @@ c
operation = 'rt_tddft'
else if (inp_compare(.false.,test,'gw')) then
operation = 'gw'
gwbse = .true.
else if (inp_compare(.false.,test,'bse')) then
operation = 'bse'
gwbse = .true.
else if (inp_compare(.false.,test,'fcidump')) then
operation = 'fcidump'
c else if (inp_compare(.false.,test,'prepare')) then
@ -377,6 +382,7 @@ c
c
if (theory .eq. ' ') then
if(osmd) theory="smd"
if(gwbse) theory='dft'
end if
if (theory .ne. ' ') then
if (.not. rtdb_cput(rtdb, 'task:theory', 1, theory))

View file

@ -769,6 +769,7 @@ $(STAMP_DIR):
clean:
@test -e $(BUILDDIR)/Makefile && { cd $(BUILDDIR) && $(MAKE) clean; } || echo "Not configured"
rm -f build install
rm -f ../mpi_include.txt ../ga*txt
realclean:

View file

@ -16,7 +16,7 @@ c (c) Battelle, PNNL, 2002.
write(LuOut,*)
write(LuOut,*) title
write(LuOut,9000)
do k = 0,mdim/6
do k = 0,(mdim-1)/6
write(LuOut,*)
write(LuOut,9010) (j,j=k*6+1,min(mdim,k*6+6))
do i=1,ndim

View file

@ -31,6 +31,7 @@ check_file () {
grep -s 'l DFT energy' $file |tail
echo ============================================================
grep -s 'Total PSPW energy' $file |tail
grep -s -B 3 "Root 1 singlet b1u" $file |grep E-0
fi
}
@ -46,6 +47,7 @@ check_file h2o2-response
check_file pspw
check_file pspw_md
check_file aump2
check_file bse_ethene
# stuff for OpenBLAS issues on macos
echo 'looking for cgemm_incopy.o'
find $TRAVIS_BUILD_DIR/src/libext -name cgemm_incopy.o

View file

@ -28,6 +28,7 @@ fi
echo NWCHEM_TOP is $NWCHEM_TOP
#TARBALL=https://github.com/nwchemgit/nwchem/releases/download/v7.0.0-beta1/nwchem-7.0.0-release.revision-5bcf0416-src.2019-11-01.tar.bz2
export USE_MPI=y
export MPICH_INTERFACE_HOSTNAME=localhost # fix for Mpich interface errors
if [[ "$FC" == "flang" ]]; then
export BUILD_MPICH=1
fi
@ -184,9 +185,9 @@ if [[ -z "$USE_INTERNALBLAS" ]]; then
# if [[ ${FC} == gfortran ]] ; then
if [[ "$MPI_IMPL" == "mpich" ]]; then
export MPICH_FC=$FC
if [[ "$arch" != "aarch64" ]]; then
export BUILD_MPICH=1
fi
# if [[ "$arch" != "aarch64" ]]; then
# export BUILD_MPICH=1
# fi
fi
if [[ `${CC} -dM -E - < /dev/null 2> /dev/null | grep -c clang` == 0 ]] && [[ `${CC} -dM -E - < /dev/null 2> /dev/null | grep -c GNU` > 0 ]] && [[ "$(expr `${CC} -dumpversion | cut -f1 -d.` \> 7)" == 1 ]]; then
if [[ "$os" == "Linux" ]] && [[ "$arch" == "x86_64" ]]; then

View file

@ -35,7 +35,7 @@ if [[ ! -z "$USE_OPENMP" ]]; then
export OMP_NUM_THREADS="$USE_OPENMP"
export OMP_STACKSIZE=32M
fi
if [[ "$os" == "Darwin" && "$BUILD_MPICH" == 1 ]] || [[ "$os" == "Darwin" && "$ARMCI_NETWORK" == MPI-PT ]] ; then
if [[ "$os" == "Darwin" && "$ARMCI_NETWORK" == MPI-PT ]] ; then
nprocs=1
fi
env|egrep MP
@ -88,6 +88,14 @@ env|egrep MP
;;
esac
;;
ARMCI)
case "$os" in
# fixes segfault for bse_ethene in ga_copy_patch...ARMCII_Iov_op_datatype
Linux)
ulimit -s unlimited
;;
esac
;;
esac
if [[ "$MPI_IMPL" == "openmpi" ]]; then
export MPIRUN_NPOPT=" --allow-run-as-root -mca mpi_yield_when_idle 0 --oversubscribe -np "
@ -132,6 +140,9 @@ fi
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs ritddft_h2o ritddft_co
cd $TRAVIS_BUILD_DIR/QA && NWCHEM_BASIS_LIBRARY=${NWCHEM_TOP}/src/basis/libraries.bse/ ./runtests.mpi.unix procs $nprocs gw_closedshell gw_openshell gw_symmetry
fi
if [[ ! $(grep -i bse $TRAVIS_BUILD_DIR/src/stubs.F| awk '/bse_input/') ]]; then
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs bse_ethene
fi
if [[ ! $(grep -i dplot $TRAVIS_BUILD_DIR/src/stubs.F| awk '/dplot_input/') ]]; then
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs rt_tddft_dplot
fi