From 030435f90c1a80867c270271f4cc1980aa3f0cfb Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 1 Jun 2020 11:35:17 +0200 Subject: [PATCH] non dTDA --- GoDuck | 1 - input/basis | 125 +++++------ input/weight | 206 ++++++++++++------ src/QuAcK/ACFDT_correlation_energy.f90 | 14 +- src/QuAcK/Bethe_Salpeter.f90 | 34 ++- .../Bethe_Salpeter_AB_matrix_dynamic.f90 | 111 ++++++++++ src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 | 76 ------- .../Bethe_Salpeter_ZAB_matrix_dynamic.f90 | 94 ++++++++ .../Bethe_Salpeter_ZB_matrix_dynamic.f90 | 66 ------ .../Bethe_Salpeter_dynamic_perturbation.f90 | 68 +++--- ...alpeter_dynamic_perturbation_iterative.f90 | 39 ++-- src/QuAcK/orthogonalization_matrix.f90 | 4 +- 12 files changed, 500 insertions(+), 338 deletions(-) create mode 100644 src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 delete mode 100644 src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 create mode 100644 src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 delete mode 100644 src/QuAcK/Bethe_Salpeter_ZB_matrix_dynamic.f90 diff --git a/GoDuck b/GoDuck index 97684a7..13eb269 100755 --- a/GoDuck +++ b/GoDuck @@ -12,7 +12,6 @@ then cp examples/molecule."$1" input/molecule cp examples/basis."$1"."$2" input/basis cp basis/"$2" input/basis.qcaml - cp examples/basis."$1"."$2" input/weight # ./bin/IntPak ./bin/QuAcK fi diff --git a/input/basis b/input/basis index d31581b..bbe0bfe 100644 --- a/input/basis +++ b/input/basis @@ -1,92 +1,71 @@ -1 14 +1 9 S 8 - 1 11420.0000000 0.0005230 - 2 1712.0000000 0.0040450 - 3 389.3000000 0.0207750 - 4 110.0000000 0.0807270 - 5 35.5700000 0.2330740 - 6 12.5400000 0.4335010 - 7 4.6440000 0.3474720 - 8 0.5118000 -0.0085080 + 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 11420.0000000 -0.0001150 - 2 1712.0000000 -0.0008950 - 3 389.3000000 -0.0046240 - 4 110.0000000 -0.0185280 - 5 35.5700000 -0.0573390 - 6 12.5400000 -0.1320760 - 7 4.6440000 -0.1725100 - 8 0.5118000 0.5999440 + 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 1.2930000 1.0000000 + 1 0.2248000 1.0000000 S 1 - 1 0.1787000 1.0000000 -S 1 - 1 0.0576000 1.0000000 + 1 0.0612400 1.0000000 P 3 - 1 26.6300000 0.0146700 - 2 5.9480000 0.0917640 - 3 1.7420000 0.2986830 + 1 13.5500000 0.0399190 + 2 2.9170000 0.2171690 + 3 0.7973000 0.5103190 P 1 - 1 0.5550000 1.0000000 + 1 0.2185000 1.0000000 P 1 - 1 0.1725000 1.0000000 -P 1 - 1 0.0491000 1.0000000 + 1 0.0561100 1.0000000 D 1 - 1 1.6540000 1.0000000 + 1 0.8170000 1.0000000 D 1 - 1 0.4690000 1.0000000 -D 1 - 1 0.1510000 1.0000000 -F 1 - 1 1.0930000 1.0000000 -F 1 - 1 0.3640000 1.0000000 -2 14 + 1 0.2300000 1.0000000 +2 9 S 8 - 1 11420.0000000 0.0005230 - 2 1712.0000000 0.0040450 - 3 389.3000000 0.0207750 - 4 110.0000000 0.0807270 - 5 35.5700000 0.2330740 - 6 12.5400000 0.4335010 - 7 4.6440000 0.3474720 - 8 0.5118000 -0.0085080 + 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 11420.0000000 -0.0001150 - 2 1712.0000000 -0.0008950 - 3 389.3000000 -0.0046240 - 4 110.0000000 -0.0185280 - 5 35.5700000 -0.0573390 - 6 12.5400000 -0.1320760 - 7 4.6440000 -0.1725100 - 8 0.5118000 0.5999440 + 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 1.2930000 1.0000000 + 1 0.2248000 1.0000000 S 1 - 1 0.1787000 1.0000000 -S 1 - 1 0.0576000 1.0000000 + 1 0.0612400 1.0000000 P 3 - 1 26.6300000 0.0146700 - 2 5.9480000 0.0917640 - 3 1.7420000 0.2986830 + 1 13.5500000 0.0399190 + 2 2.9170000 0.2171690 + 3 0.7973000 0.5103190 P 1 - 1 0.5550000 1.0000000 + 1 0.2185000 1.0000000 P 1 - 1 0.1725000 1.0000000 -P 1 - 1 0.0491000 1.0000000 + 1 0.0561100 1.0000000 D 1 - 1 1.6540000 1.0000000 + 1 0.8170000 1.0000000 D 1 - 1 0.4690000 1.0000000 -D 1 - 1 0.1510000 1.0000000 -F 1 - 1 1.0930000 1.0000000 -F 1 - 1 0.3640000 1.0000000 - + 1 0.2300000 1.0000000 diff --git a/input/weight b/input/weight index d31581b..5dc5a6c 100644 --- a/input/weight +++ b/input/weight @@ -1,92 +1,174 @@ 1 14 S 8 - 1 11420.0000000 0.0005230 - 2 1712.0000000 0.0040450 - 3 389.3000000 0.0207750 - 4 110.0000000 0.0807270 - 5 35.5700000 0.2330740 - 6 12.5400000 0.4335010 - 7 4.6440000 0.3474720 - 8 0.5118000 -0.0085080 + 1 8236.0000000 0.0005310 + 2 1235.0000000 0.0041080 + 3 280.8000000 0.0210870 + 4 79.2700000 0.0818530 + 5 25.5900000 0.2348170 + 6 8.9970000 0.4344010 + 7 3.3190000 0.3461290 + 8 0.3643000 -0.0089830 S 8 - 1 11420.0000000 -0.0001150 - 2 1712.0000000 -0.0008950 - 3 389.3000000 -0.0046240 - 4 110.0000000 -0.0185280 - 5 35.5700000 -0.0573390 - 6 12.5400000 -0.1320760 - 7 4.6440000 -0.1725100 - 8 0.5118000 0.5999440 + 1 8236.0000000 -0.0001130 + 2 1235.0000000 -0.0008780 + 3 280.8000000 -0.0045400 + 4 79.2700000 -0.0181330 + 5 25.5900000 -0.0557600 + 6 8.9970000 -0.1268950 + 7 3.3190000 -0.1703520 + 8 0.3643000 0.5986840 S 1 - 1 1.2930000 1.0000000 + 1 0.9059000 1.0000000 S 1 - 1 0.1787000 1.0000000 + 1 0.1285000 1.0000000 S 1 - 1 0.0576000 1.0000000 + 1 0.0440200 1.0000000 P 3 - 1 26.6300000 0.0146700 - 2 5.9480000 0.0917640 - 3 1.7420000 0.2986830 + 1 18.7100000 0.0140310 + 2 4.1330000 0.0868660 + 3 1.2000000 0.2902160 P 1 - 1 0.5550000 1.0000000 + 1 0.3827000 1.0000000 P 1 - 1 0.1725000 1.0000000 + 1 0.1209000 1.0000000 P 1 - 1 0.0491000 1.0000000 + 1 0.0356900 1.0000000 D 1 - 1 1.6540000 1.0000000 + 1 1.0970000 1.0000000 D 1 - 1 0.4690000 1.0000000 + 1 0.3180000 1.0000000 D 1 - 1 0.1510000 1.0000000 + 1 0.1000000 1.0000000 F 1 - 1 1.0930000 1.0000000 + 1 0.7610000 1.0000000 F 1 - 1 0.3640000 1.0000000 + 1 0.2680000 1.0000000 2 14 S 8 - 1 11420.0000000 0.0005230 - 2 1712.0000000 0.0040450 - 3 389.3000000 0.0207750 - 4 110.0000000 0.0807270 - 5 35.5700000 0.2330740 - 6 12.5400000 0.4335010 - 7 4.6440000 0.3474720 - 8 0.5118000 -0.0085080 + 1 8236.0000000 0.0005310 + 2 1235.0000000 0.0041080 + 3 280.8000000 0.0210870 + 4 79.2700000 0.0818530 + 5 25.5900000 0.2348170 + 6 8.9970000 0.4344010 + 7 3.3190000 0.3461290 + 8 0.3643000 -0.0089830 S 8 - 1 11420.0000000 -0.0001150 - 2 1712.0000000 -0.0008950 - 3 389.3000000 -0.0046240 - 4 110.0000000 -0.0185280 - 5 35.5700000 -0.0573390 - 6 12.5400000 -0.1320760 - 7 4.6440000 -0.1725100 - 8 0.5118000 0.5999440 + 1 8236.0000000 -0.0001130 + 2 1235.0000000 -0.0008780 + 3 280.8000000 -0.0045400 + 4 79.2700000 -0.0181330 + 5 25.5900000 -0.0557600 + 6 8.9970000 -0.1268950 + 7 3.3190000 -0.1703520 + 8 0.3643000 0.5986840 S 1 - 1 1.2930000 1.0000000 + 1 0.9059000 1.0000000 S 1 - 1 0.1787000 1.0000000 + 1 0.1285000 1.0000000 S 1 - 1 0.0576000 1.0000000 + 1 0.0440200 1.0000000 P 3 - 1 26.6300000 0.0146700 - 2 5.9480000 0.0917640 - 3 1.7420000 0.2986830 + 1 18.7100000 0.0140310 + 2 4.1330000 0.0868660 + 3 1.2000000 0.2902160 P 1 - 1 0.5550000 1.0000000 + 1 0.3827000 1.0000000 P 1 - 1 0.1725000 1.0000000 + 1 0.1209000 1.0000000 P 1 - 1 0.0491000 1.0000000 + 1 0.0356900 1.0000000 D 1 - 1 1.6540000 1.0000000 + 1 1.0970000 1.0000000 D 1 - 1 0.4690000 1.0000000 + 1 0.3180000 1.0000000 D 1 - 1 0.1510000 1.0000000 + 1 0.1000000 1.0000000 F 1 - 1 1.0930000 1.0000000 + 1 0.7610000 1.0000000 F 1 - 1 0.3640000 1.0000000 - - + 1 0.2680000 1.0000000 +3 9 +S 3 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 +S 1 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 +P 1 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +4 9 +S 3 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 +S 1 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 +P 1 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +5 9 +S 3 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 +S 1 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 +P 1 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +6 9 +S 3 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 +S 1 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 +P 1 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 diff --git a/src/QuAcK/ACFDT_correlation_energy.f90 b/src/QuAcK/ACFDT_correlation_energy.f90 index 4dead7b..82b578a 100644 --- a/src/QuAcK/ACFDT_correlation_energy.f90 +++ b/src/QuAcK/ACFDT_correlation_energy.f90 @@ -69,8 +69,18 @@ subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ER ! Compute Tr(K x P_lambda) - X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:)) - Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:)) +! X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:)) +! Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:)) + +! print*,'X' +! call matout(nS,nS,X) +! print*,'Y' +! call matout(nS,nS,Y) + +! print*,'Ap' +! call matout(nS,nS,Ap) +! print*,'Bp' +! call matout(nS,nS,Bp) EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) & + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index b73504d..9d19585 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -26,6 +26,7 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & ! Local variables logical :: evDyn = .false. + logical :: W_BSE = .false. integer :: ispin double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: XpY_BSE(:,:,:) @@ -39,7 +40,8 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & ! Memory allocation - allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin),rho_BSE(nBas,nBas,nS,nspin)) + allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + if(W_BSE) allocate(rho_BSE(nBas,nBas,nS,nspin)) !------------------- ! Singlet manifold @@ -63,7 +65,10 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) -! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + + ! Compute the dynamical screening at the BSE level + + if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -73,8 +78,14 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + if(W_BSE) then + call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + else + call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + end if + end if end if @@ -101,7 +112,10 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) -! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + + ! Compute the dynamical screening at the BSE level + + if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -111,8 +125,14 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + if(W_BSE) then + call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + else + call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + end if + end if end if diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 new file mode 100644 index 0000000..d41f252 --- /dev/null +++ b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -0,0 +1,111 @@ +subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,Ap,Am,Bp,Bm) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: OmBSE + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: maxS + double precision :: chi_A,chi_B,eps + double precision :: chi_Ap,chi_Am,chi_Bp,chi_Bm + double precision :: eps_Ap,eps_Am,eps_Bp,eps_Bm + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: Ap(nS,nS) + double precision,intent(out) :: Am(nS,nS) + double precision,intent(out) :: Bp(nS,nS) + double precision,intent(out) :: Bm(nS,nS) + +! Initialization + + Ap(:,:) = 0d0 + Am(:,:) = 0d0 + Bp(:,:) = 0d0 + Bm(:,:) = 0d0 + +! Number of poles taken into account + + maxS = nS + +! Build dynamic A matrix + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi_A = 0d0 + chi_B = 0d0 + + do kc=1,maxS + + eps = OmRPA(kc)**2 + eta**2 + chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/eps + chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/eps + + enddo + + Ap(ia,jb) = Ap(ia,jb) - 4d0*lambda*chi_A + Am(ia,jb) = Am(ia,jb) - 4d0*lambda*chi_A + Bp(ia,jb) = Bp(ia,jb) - 4d0*lambda*chi_B + Bm(ia,jb) = Bm(ia,jb) - 4d0*lambda*chi_B + + chi_Ap = 0d0 + chi_Am = 0d0 + chi_Bp = 0d0 + chi_Bm = 0d0 + + do kc=1,maxS + + eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 + eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am + + eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 + eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am + + eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 + eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm + + eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 + eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm + + enddo + + Ap(ia,jb) = Ap(ia,jb) - 2d0*lambda*chi_Ap + Am(ia,jb) = Am(ia,jb) - 2d0*lambda*chi_Am + + Bp(ia,jb) = Bp(ia,jb) - 2d0*lambda*chi_Bp + Bm(ia,jb) = Bm(ia,jb) - 2d0*lambda*chi_Bm + + enddo + enddo + enddo + enddo + +end subroutine Bethe_Salpeter_AB_matrix_dynamic diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 deleted file mode 100644 index 3e93634..0000000 --- a/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 +++ /dev/null @@ -1,76 +0,0 @@ -subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,B_dyn) - -! Compute the dynamic part of the Bethe-Salpeter equation matrices - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: eta - double precision,intent(in) :: lambda - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: OmRPA(nS) - double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: maxS - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: B_dyn(nS,nS) - -! Initialization - - B_dyn(:,:) = 0d0 - -! Number of poles taken into account - - maxS = nS - -! Build dynamic A matrix - - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - - chi = 0d0 - do kc=1,maxS - - eps = OmRPA(kc)**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/eps - - enddo - - B_dyn(ia,jb) = B_dyn(ia,jb) - 4d0*lambda*chi - - chi = 0d0 - do kc=1,maxS - - eps = (OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps - - eps = (OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps - - enddo - - B_dyn(ia,jb) = B_dyn(ia,jb) - 2d0*lambda*chi - - enddo - enddo - enddo - enddo - -end subroutine Bethe_Salpeter_B_matrix_dynamic diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 new file mode 100644 index 0000000..da0db60 --- /dev/null +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -0,0 +1,94 @@ +subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZAp,ZAm,ZBp,ZBm) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: OmBSE + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: maxS + double precision :: chi_Ap,chi_Am,chi_Bp,chi_Bm + double precision :: eps_Ap,eps_Am,eps_Bp,eps_Bm + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: ZAp(nS,nS) + double precision,intent(out) :: ZAm(nS,nS) + double precision,intent(out) :: ZBp(nS,nS) + double precision,intent(out) :: ZBm(nS,nS) + +! Initialization + + ZAp(:,:) = 0d0 + ZAm(:,:) = 0d0 + ZBp(:,:) = 0d0 + ZBm(:,:) = 0d0 + +! Number of poles taken into account + + maxS = nS + +! Build dynamic A matrix + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi_Ap = 0d0 + chi_Am = 0d0 + chi_Bp = 0d0 + chi_Bm = 0d0 + + do kc=1,maxS + + eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 + eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap)**2 + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am)**2 + + eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 + eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 + chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap)**2 + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am)**2 + + eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 + eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp)**2 + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm)**2 + + eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 + eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 + chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp)**2 + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm)**2 + + enddo + + ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap + ZAm(ia,jb) = ZAm(ia,jb) + 2d0*lambda*chi_Am + + ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp + ZBm(ia,jb) = ZBm(ia,jb) + 2d0*lambda*chi_Bm + + enddo + enddo + enddo + enddo + +end subroutine Bethe_Salpeter_ZAB_matrix_dynamic diff --git a/src/QuAcK/Bethe_Salpeter_ZB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZB_matrix_dynamic.f90 deleted file mode 100644 index a337b5e..0000000 --- a/src/QuAcK/Bethe_Salpeter_ZB_matrix_dynamic.f90 +++ /dev/null @@ -1,66 +0,0 @@ -subroutine Bethe_Salpeter_ZB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZB_dyn) - -! Compute the dynamic part of the Bethe-Salpeter equation matrices - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: eta - double precision,intent(in) :: lambda - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: OmRPA(nS) - double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: maxS - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: ZB_dyn(nS,nS) - -! Initialization - - ZB_dyn(:,:) = 0d0 - -! Number of poles taken into account - - maxS = nS - -! Build dynamic A matrix - - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - - chi = 0d0 - do kc=1,maxS - - eps = (OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*((OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps)**2 - - eps = (OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*((OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps)**2 - - enddo - - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + 2d0*lambda*chi - - enddo - enddo - enddo - enddo - -end subroutine Bethe_Salpeter_ZB_matrix_dynamic diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index b8b396a..3e41fae 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -25,7 +25,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O ! Local variables - logical :: TDA_dyn = .true. + logical :: dTDA = .false. integer :: ia integer,parameter :: maxS = 10 double precision :: gapGW @@ -34,16 +34,22 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O double precision,allocatable :: ZDyn(:) double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: A_dyn(:,:) - double precision,allocatable :: B_dyn(:,:) - double precision,allocatable :: ZA_dyn(:,:) - double precision,allocatable :: ZB_dyn(:,:) + + double precision,allocatable :: Ap_dyn(:,:) + double precision,allocatable :: Am_dyn(:,:) + double precision,allocatable :: ZAp_dyn(:,:) + double precision,allocatable :: ZAm_dyn(:,:) + + double precision,allocatable :: Bp_dyn(:,:) + double precision,allocatable :: Bm_dyn(:,:) + double precision,allocatable :: ZBp_dyn(:,:) + double precision,allocatable :: ZBm_dyn(:,:) ! Memory allocation - allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),A_dyn(nS,nS),ZA_dyn(nS,nS)) + allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - if(.not.TDA_dyn) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS)) + if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS),ZBp_dyn(nS,nS),ZBm_dyn(nS,nS)) gapGW = eGW(nO+1) - eGW(nO) @@ -58,40 +64,42 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - ! Resonant part of the BSE correction - - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),A_dyn(:,:)) - - ! Renormalization factor of the resonant part - - call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),ZA_dyn(:,:)) - ! First-order correction - if(TDA_dyn) then + if(dTDA) then - ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) - OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) + ! Resonant part of the BSE correction for dynamical TDA + + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),Ap_dyn(:,:)) + + ! Renormalization factor of the resonant parts for dynamical TDA + + call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),ZAp_dyn(:,:)) + + ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) else - ! Anti-resonant part of the BSE correction + ! Resonant and anti-resonant part of the BSE correction - call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),B_dyn(:,:)) + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & + Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) - ! Renormalization factor of the anti-resonant part + ! Renormalization factor of the resonant and anti-resonant parts - call Bethe_Salpeter_ZB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),ZB_dyn(:,:)) + call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & + ZAp_dyn(:,:),ZAm_dyn(:,:),ZBp_dyn(:,:),ZBm_dyn(:,:)) - ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(ZA_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(ZB_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(ZB_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:))) - OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(B_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) end if diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index 35777c1..2cd4eea 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 @@ -25,7 +25,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n ! Local variables - logical :: TDA_dyn = .true. + logical :: dTDA = .true. integer :: ia integer,parameter :: maxS = 10 double precision :: gapGW @@ -40,16 +40,16 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n double precision,allocatable :: OmOld(:) double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: A_dyn(:,:) - double precision,allocatable :: B_dyn(:,:) - double precision,allocatable :: ZA_dyn(:,:) - double precision,allocatable :: ZB_dyn(:,:) + double precision,allocatable :: Ap_dyn(:,:) + double precision,allocatable :: Am_dyn(:,:) + double precision,allocatable :: Bp_dyn(:,:) + double precision,allocatable :: Bm_dyn(:,:) ! Memory allocation - allocate(OmDyn(nS),OmOld(nS),X(nS),Y(nS),A_dyn(nS,nS)) + allocate(OmDyn(nS),OmOld(nS),X(nS),Y(nS),Ap_dyn(nS,nS)) - if(.not.TDA_dyn) allocate(B_dyn(nS,nS)) + if(.not.dTDA) allocate(Am_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS)) gapGW = eGW(nO+1) - eGW(nO) @@ -78,27 +78,28 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - ! Resonant part of the BSE correction - - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:),A_dyn(:,:)) - ! First-order correction - if(TDA_dyn) then + if(dTDA) then + + ! Resonant part of the BSE correction - OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:),Ap_dyn(:,:)) + + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) else ! Anti-resonant part of the BSE correction - call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:),B_dyn(:,:)) - - OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(B_dyn(:,:),X(:))) + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & + Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) + end if write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6)') & diff --git a/src/QuAcK/orthogonalization_matrix.f90 b/src/QuAcK/orthogonalization_matrix.f90 index e60f253..348abdb 100644 --- a/src/QuAcK/orthogonalization_matrix.f90 +++ b/src/QuAcK/orthogonalization_matrix.f90 @@ -13,7 +13,7 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X) logical :: debug double precision,allocatable :: UVec(:,:),Uval(:) - double precision,parameter :: thresh = 1d-8 + double precision,parameter :: thresh = 1d-6 integer :: i @@ -45,7 +45,7 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X) if(Uval(i) < thresh) then - write(*,*) 'Eigenvalue',i,'is very small in Lowdin orthogonalization = ',Uval(i) + write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i) endif