PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-03 18:16:03 +01:00
Code Releases Activity

(1)

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-06-05 22:34:32 +02:00
parent bf7f441d4d
commit 769f4e9221
5 changed files with 104 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

127
input/basis
View File

@ -1,71 +1,64 @@
1 9 1 6
S 8 S 9
1 9046.0000000 0.0007000 1 9.046000E+03 7.000000E-04
2 1357.0000000 0.0053890 2 1.357000E+03 5.389000E-03
3 309.3000000 0.0274060 3 3.093000E+02 2.740600E-02
4 87.7300000 0.1032070 4 8.773000E+01 1.032070E-01
5 28.5600000 0.2787230 5 2.856000E+01 2.787230E-01
6 10.2100000 0.4485400 6 1.021000E+01 4.485400E-01
7 3.8380000 0.2782380 7 3.838000E+00 2.782380E-01
8 0.7466000 0.0154400 8 7.466000E-01 1.544000E-02
S 8 9 2.248000E-01 -2.864000E-03
1 9046.0000000 -0.0001530 S 9
2 1357.0000000 -0.0012080 1 9.046000E+03 -1.530000E-04
3 309.3000000 -0.0059920 2 1.357000E+03 -1.208000E-03
4 87.7300000 -0.0245440 3 3.093000E+02 -5.992000E-03
5 28.5600000 -0.0674590 4 8.773000E+01 -2.454400E-02
6 10.2100000 -0.1580780 5 2.856000E+01 -6.745900E-02
7 3.8380000 -0.1218310 6 1.021000E+01 -1.580780E-01
8 0.7466000 0.5490030 7 3.838000E+00 -1.218310E-01
8 7.466000E-01 5.490030E-01
9 2.248000E-01 5.788150E-01
S 1 S 1
1 0.2248000 1.0000000 1 2.248000E-01 1.000000E+00
P 4
1 1.355000E+01 3.991900E-02
2 2.917000E+00 2.171690E-01
3 7.973000E-01 5.103190E-01
4 2.185000E-01 4.622140E-01
P 1
1 2.185000E-01 1.000000E+00
D 1
1 8.170000E-01 1.0000000
2 6
S 9
1 9.046000E+03 7.000000E-04
2 1.357000E+03 5.389000E-03
3 3.093000E+02 2.740600E-02
4 8.773000E+01 1.032070E-01
5 2.856000E+01 2.787230E-01
6 1.021000E+01 4.485400E-01
7 3.838000E+00 2.782380E-01
8 7.466000E-01 1.544000E-02
9 2.248000E-01 -2.864000E-03
S 9
1 9.046000E+03 -1.530000E-04
2 1.357000E+03 -1.208000E-03
3 3.093000E+02 -5.992000E-03
4 8.773000E+01 -2.454400E-02
5 2.856000E+01 -6.745900E-02
6 1.021000E+01 -1.580780E-01
7 3.838000E+00 -1.218310E-01
8 7.466000E-01 5.490030E-01
9 2.248000E-01 5.788150E-01
S 1 S 1
1 0.0612400 1.0000000 1 2.248000E-01 1.000000E+00
P 3 P 4
1 13.5500000 0.0399190 1 1.355000E+01 3.991900E-02
2 2.9170000 0.2171690 2 2.917000E+00 2.171690E-01
3 0.7973000 0.5103190 3 7.973000E-01 5.103190E-01
4 2.185000E-01 4.622140E-01
P 1 P 1
1 0.2185000 1.0000000 1 2.185000E-01 1.000000E+00
P 1
1 0.0561100 1.0000000
D 1 D 1
1 0.8170000 1.0000000 1 8.170000E-01 1.0000000
D 1
1 0.2300000 1.0000000
2 9
S 8
1 9046.0000000 0.0007000
2 1357.0000000 0.0053890
3 309.3000000 0.0274060
4 87.7300000 0.1032070
5 28.5600000 0.2787230
6 10.2100000 0.4485400
7 3.8380000 0.2782380
8 0.7466000 0.0154400
S 8
1 9046.0000000 -0.0001530
2 1357.0000000 -0.0012080
3 309.3000000 -0.0059920
4 87.7300000 -0.0245440
5 28.5600000 -0.0674590
6 10.2100000 -0.1580780
7 3.8380000 -0.1218310
8 0.7466000 0.5490030
S 1
1 0.2248000 1.0000000
S 1
1 0.0612400 1.0000000
P 3
1 13.5500000 0.0399190
2 2.9170000 0.2171690
3 0.7973000 0.5103190
P 1
1 0.2185000 1.0000000
P 1
1 0.0561100 1.0000000
D 1
1 0.8170000 1.0000000
D 1
1 0.2300000 1.0000000

4
input/methods
View File

@ -9,11 +9,11 @@
# CIS CID CISD # CIS CID CISD
F F F F F F
# RPA RPAx ppRPA # RPA RPAx ppRPA
F F F F T F
# G0F2 evGF2 G0F3 evGF3 # G0F2 evGF2 G0F3 evGF3
F F F F F F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW
T F F F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F F F F
# MCMP2 # MCMP2

38
src/QuAcK/QuAcK.f90
View File

@ -383,7 +383,7 @@ program QuAcK
if(doCCD) then if(doCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CCD)
@ -402,7 +402,7 @@ program QuAcK
if(doCCSD) then if(doCCSD) then
call cpu_time(start_CCSD) call cpu_time(start_CCSD)
call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC(1),nO(1),nV(1),nR(1), & call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCSD) call cpu_time(end_CCSD)
@ -419,7 +419,7 @@ program QuAcK
if(do_drCCD) then if(do_drCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CCD)
@ -436,7 +436,7 @@ program QuAcK
if(do_rCCD) then if(do_rCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CCD)
@ -453,7 +453,7 @@ program QuAcK
if(do_lCCD) then if(do_lCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CCD)
@ -470,7 +470,7 @@ program QuAcK
if(do_pCCD) then if(do_pCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CCD)
@ -570,7 +570,7 @@ program QuAcK
call cpu_time(start_ppRPA) call cpu_time(start_ppRPA)
call ppRPA(singlet_manifold,triplet_manifold, & call ppRPA(singlet_manifold,triplet_manifold, &
nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO,eHF) nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA) call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA t_ppRPA = end_ppRPA - start_ppRPA
@ -587,7 +587,7 @@ program QuAcK
! call cpu_time(start_ADC) ! call cpu_time(start_ADC)
! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & ! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, &
! nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO) ! nBas,nC,nO,nV,nR,eHF,ERI_MO)
! call cpu_time(end_ADC) ! call cpu_time(end_ADC)
! t_ADC = end_ADC - start_ADC ! t_ADC = end_ADC - start_ADC
@ -604,7 +604,7 @@ program QuAcK
call cpu_time(start_GF2) call cpu_time(start_GF2)
call G0F2(BSE_GF,TDA_GF,singlet_manifold,triplet_manifold,linGF, & call G0F2(BSE_GF,TDA_GF,singlet_manifold,triplet_manifold,linGF, &
eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_GF2) call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2 t_GF2 = end_GF2 - start_GF2
@ -621,7 +621,7 @@ program QuAcK
call cpu_time(start_GF2) call cpu_time(start_GF2)
call evGF2(BSE_GF,TDA_GF,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, & call evGF2(BSE_GF,TDA_GF,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, &
eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_GF2) call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2 t_GF2 = end_GF2 - start_GF2
@ -637,7 +637,7 @@ program QuAcK
if(doG0F3) then if(doG0F3) then
call cpu_time(start_GF3) call cpu_time(start_GF3)
call G0F3(renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) call G0F3(renormGF,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_GF3) call cpu_time(end_GF3)
t_GF3 = end_GF3 - start_GF3 t_GF3 = end_GF3 - start_GF3
@ -653,7 +653,7 @@ program QuAcK
if(doevGF3) then if(doevGF3) then
call cpu_time(start_GF3) call cpu_time(start_GF3)
call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_GF3) call cpu_time(end_GF3)
t_GF3 = end_GF3 - start_GF3 t_GF3 = end_GF3 - start_GF3
@ -673,7 +673,7 @@ program QuAcK
call cpu_time(start_G0W0) call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, & call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0) call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0 t_G0W0 = end_G0W0 - start_G0W0
@ -691,7 +691,7 @@ program QuAcK
call cpu_time(start_evGW) call cpu_time(start_evGW)
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_evGW) call cpu_time(end_evGW)
t_evGW = end_evGW - start_evGW t_evGW = end_evGW - start_evGW
@ -709,7 +709,7 @@ program QuAcK
call cpu_time(start_qsGW) call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF)
call cpu_time(end_qsGW) call cpu_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW t_qsGW = end_qsGW - start_qsGW
@ -729,7 +729,7 @@ program QuAcK
call cpu_time(start_G0T0) call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, & call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_G0T0) call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0 t_G0T0 = end_G0T0 - start_G0T0
@ -747,7 +747,7 @@ program QuAcK
call cpu_time(start_evGT) call cpu_time(start_evGT)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
BSE_GW,TDA_GW,singlet_manifold,triplet_manifold,eta_GW, & BSE_GW,TDA_GW,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_evGT) call cpu_time(end_evGT)
t_evGT = end_evGT - start_evGT t_evGT = end_evGT - start_evGT
@ -864,7 +864,7 @@ program QuAcK
call cpu_time(start_G0W0) call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, & call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0) call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0 t_G0W0 = end_G0W0 - start_G0W0
@ -878,7 +878,7 @@ program QuAcK
call cpu_time(start_G0T0) call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, & call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0)
call cpu_time(end_G0T0) call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0 t_G0T0 = end_G0T0 - start_G0T0

3
src/QuAcK/RPAx.f90
View File

@ -62,6 +62,8 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
@ -74,6 +76,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif

32
src/QuAcK/print_transition_vectors.f90
View File

@ -22,6 +22,7 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
integer :: ia,jb,i,j,a,b integer :: ia,jb,i,j,a,b
integer,parameter :: maxS = 10 integer,parameter :: maxS = 10
double precision,parameter :: thres_vec = 0.1d0 double precision,parameter :: thres_vec = 0.1d0
double precision :: norm
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
@ -32,31 +33,38 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
write(*,*) write(*,*)
do ia=1,min(nS,maxS) do ia=1,min(nS,maxS)
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
norm = 0d0
do jb=1,nS
norm = norm + X(jb)*X(jb)
end do
norm = sqrt(norm)
print*,'--------------------------------' print*,'--------------------------------'
write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV' write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV'
print*,'--------------------------------' print*,'--------------------------------'
print*,'Resonant vectors'
print*,'--------------------------------'
jb = 0 jb = 0
do j=nC+1,nO do j=nC+1,nO
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',X(jb) if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0)
end do end do
end do end do
print*,'--------------------------------' norm = 0d0
print*,'Anti-resonant vectors' do jb=1,nS
print*,'--------------------------------' norm = norm + Y(jb)*Y(jb)
end do
norm = sqrt(norm)
jb = 0 jb = 0
do j=nC+1,nO do j=nC+1,nO
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',Y(jb) if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0)
end do end do
end do end do
write(*,*) write(*,*)