diff --git a/src/prg_ewald_mod.F90 b/src/prg_ewald_mod.F90 index 5539c12d..a56de96e 100644 --- a/src/prg_ewald_mod.F90 +++ b/src/prg_ewald_mod.F90 @@ -436,8 +436,8 @@ subroutine Ewald_Real_Space_Test(COULOMBV,I,RX,RY,RZ,LBox, & Ra(2) = RY(I) Ra(3) = RZ(I) -! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nnI,J,Rb,Rab,dR,MAGR,MAGR2,TJ,DC,Z,NUMREP_ERFC,CA) & -! !$OMP REDUCTION(+:COULOMBV) + ! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nnI,J,Rb,Rab,dR,MAGR,MAGR2,TJ,DC,Z,NUMREP_ERFC,CA) & + ! !$OMP REDUCTION(+:COULOMBV) do nnI = 1,nrnnlist(I) Rb(1) = nnRx(I,nnI) Rb(2) = nnRy(I,nnI) @@ -488,7 +488,7 @@ subroutine Ewald_Real_Space_Test(COULOMBV,I,RX,RY,RZ,LBox, & endif enddo -! !$OMP END PARALLEL DO + ! !$OMP END PARALLEL DO COULOMBV = KECONST*COULOMBV end subroutine Ewald_Real_Space_Test @@ -815,8 +815,8 @@ subroutine Ewald_k_Space_latte_single(COULOMBV,J,RXYZ,Box,DELTAQ,Nr_atoms,COULAC M22 = M*RECIPVECS(2,2) do N = NMIN,NMAX - K(1) = L11 - K(2) = M22 + K(1) = L11 + K(2) = M22 K(3) = N*RECIPVECS(3,3) K2 = K(1)*K(1) + K(2)*K(2) + K(3)*K(3) @@ -824,7 +824,7 @@ subroutine Ewald_k_Space_latte_single(COULOMBV,J,RXYZ,Box,DELTAQ,Nr_atoms,COULAC KEPREF = KECONST*PREFACTOR JDOT = K(1)*RXYZ(1,J) + K(2)*RXYZ(2,J) + K(3)*RXYZ(3,J) SINJDOT = sin(JDOT) - COSJDOT = cos(JDOT) + COSJDOT = cos(JDOT) do I = 1,Nr_atoms IDOT = K(1)*RXYZ(1,I) + K(2)*RXYZ(2,I) + K(3)*RXYZ(3,I) COULOMBV(I) = COULOMBV(I) + KEPREF*DELTAQ(J)*(COSJDOT*cos(IDOT)+SINJDOT*sin(IDOT)) @@ -841,9 +841,9 @@ subroutine Ewald_k_Space_latte_single(COULOMBV,J,RXYZ,Box,DELTAQ,Nr_atoms,COULAC end subroutine Ewald_k_Space_latte_single subroutine Ewald_k_Space_Test(COULOMBV,RX,RY,RZ,LBox,DELTAQ,Nr_atoms,COULACC,Max_Nr_Neigh) -! + ! implicit none -! + ! integer, parameter :: PREC = 8 real(PREC), parameter :: ONE = 1.D0, TWO = 2.D0, ZERO = 0.D0, SIX = 6.D0, THREE = 3.D0, FOUR = 4.D0 real(PREC), parameter :: FOURTYEIGHT = 48.D0, ELEVEN = 11.D0, SIXTEEN = 16.D0, EIGHT = 8.D0 @@ -864,7 +864,7 @@ subroutine Ewald_k_Space_Test(COULOMBV,RX,RY,RZ,LBox,DELTAQ,Nr_atoms,COULACC,Max real(PREC) :: PREVIR, COSSUM,SINSUM,DOT, KEPREF, COSSUM2, SINSUM2 integer :: I,J,L,M,N, ccnt, nnI, LMAX,MMAX,NMAX,NMIN,MMIN -! + ! COULVOL = LBox(1)*LBox(2)*LBox(3) SQRTX = sqrt(-log(COULACC)) @@ -946,7 +946,7 @@ subroutine Ewald_k_Space_Test(COULOMBV,RX,RY,RZ,LBox,DELTAQ,Nr_atoms,COULACC,Max SINSUM2 = SINSUM*SINSUM ! Add up energy and force contributions -! + ! KEPREF = KECONST*PREFACTOR !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I) do I = 1,Nr_atoms diff --git a/src/prg_implicit_fermi_mod.F90 b/src/prg_implicit_fermi_mod.F90 index 37d03618..20d20f7b 100644 --- a/src/prg_implicit_fermi_mod.F90 +++ b/src/prg_implicit_fermi_mod.F90 @@ -100,7 +100,7 @@ subroutine prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, call prg_conjgrad(y_bml, ai_bml, I_bml, aux_bml, d_bml, w_bml, 0.0001_dp, threshold) do i = 1, nsteps call bml_copy(I_bml, Inv_bml(i)) - enddo + enddo else ! Otherwise use previous inverse as starting guess call bml_copy(Inv_bml(1),ai_bml) @@ -123,7 +123,7 @@ subroutine prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, call bml_scale_add_identity(y_bml, 2.0_dp, 1.0_dp, threshold) ! Find inverse ai = (2*(P2-P)+I)^-1 !call prg_conjgrad(y_bml, Inv_bml(i), I_bml, aux_bml, d_bml, w_bml, tol, threshold) - !call bml_copy(Inv_bml(i),ai_bml) + !call bml_copy(Inv_bml(i),ai_bml) call prg_newtonschulz(y_bml, Inv_bml(i), d_bml, w_bml, aux_bml, I_bml, tol, threshold) call bml_multiply(Inv_bml(i), p2_bml, p_bml, 1.0_dp, 0.0_dp, threshold) !call bml_copy(ai_bml, Inv_bml(i)) ! Save inverses for use in perturbation response calculation @@ -133,37 +133,37 @@ subroutine prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, trdPdmu = beta*(trP0 - bml_sum_squares(p_bml)) ! sum p(i,j)**2 occErr = abs(trP0 - nocc) write(*,*) 'occerr =', occErr - + ! If occupation error is too large, do bisection method if (occerr > 1.0_dp) then - ! if (newerr > occerr) then + ! if (newerr > occerr) then if (nocc-trP0 < 0.0_dp .and. prev .eq. -1) then prev = -1 - else if (nocc-trP0 > 0.0_dp .and. prev .eq. 1) then + else if (nocc-trP0 > 0.0_dp .and. prev .eq. 1) then prev = 1 else if (nocc-trP0 > 0.0_dp .and. prev .eq. -1) then prev = 1 alpha = alpha/2 - ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) - ! do while(newerr > occerr) - ! alpha = alpha/2 - ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) - ! enddo - else - prev = -1 - alpha = alpha/2 - ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) - ! do while(newerr > occerr) - ! alpha = alpha/2 - ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) - ! enddo + ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) + ! do while(newerr > occerr) + ! alpha = alpha/2 + ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) + ! enddo + else + prev = -1 + alpha = alpha/2 + ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) + ! do while(newerr > occerr) + ! alpha = alpha/2 + ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) + ! enddo endif !newerr = abs(trP0+prev*alpha*trdPdmu-nocc) !do while(newerr > occerr .and. abs(occerr-newerr)/occerr<0.1 ) ! alpha = alpha/2 ! newerr = abs(trP0+prev*alpha*trdPdmu-nocc) !enddo -! write(*,*) 'newerr =', newerr + ! write(*,*) 'newerr =', newerr mu = mu + prev*alpha muadj = 1 @@ -180,9 +180,9 @@ subroutine prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, end if enddo - if (iter .ge. maxiter) then - write(*,*) 'Could not converge chemical potential in prg_impplicit_fermi_save_inverse' - end if + if (iter .ge. maxiter) then + write(*,*) 'Could not converge chemical potential in prg_impplicit_fermi_save_inverse' + end if ! Adjusting the occupation sometimes causes the perturbation calculation to not converge. ! For now we recompute the DM one extra time if mu was adjusted. !if (muadj .eq. 1) then @@ -317,9 +317,9 @@ subroutine prg_implicit_fermi(h_bml, p_bml, nsteps, k, nocc, & trdPdmu = trdPdmu - bml_sum_squares(p_bml) ! sum p(i,j)**2 trdPdmu = beta * trdPdmu occErr = abs(trP0 - nocc) - if (occErr .gt. occErrLimit) then + if (occErr .gt. occErrLimit) then mu = mu + (nocc - trP0)/trdPdmu - end if + end if write(*,*) "mu =", mu enddo @@ -503,25 +503,25 @@ subroutine prg_implicit_fermi_first_order_response(H0_bml, H1_bml, P0_bml, P1_bm call bml_multiply(Inv_bml(i), C_bml, P1_bml, 1.0_dp, 0.0_dp, threshold) enddo -! do i = 1, nsteps-1 - ! D = A^-1*P0 -! call bml_multiply(Inv_bml(i), B0_bml, C0_bml, 1.0_dp, 0.0_dp, threshold) -! call bml_multiply(C0_bml, B0_bml, B_bml, 1.0_dp, 0.0_dp, threshold) - ! B0 = A^-1*P0^2 -! call bml_copy(B_bml,B0_bml) - ! B = I + D -P0*D -! call bml_add(B_bml, C0_bml, -1.0_dp, 1.0_dp, threshold) -! call bml_scale_add_identity(B_bml, 1.0_dp, 1.0_dp, threshold) - ! P1 = 2D*P1(I+D-P0*D) -! call bml_multiply(C0_bml, P1_bml, C_bml, 1.0_dp, 0.0_dp, threshold) -! call bml_multiply(C_bml, B_bml, P1_bml, 2.0_dp, 0.0_dp, threshold) -! enddo -! call bml_multiply(B0_bml, P1_bml, C_bml, 2.0_dp, 0.0_dp, threshold) -! call bml_copy(P1_bml, B_bml) -! call bml_add(B_bml, C_bml, 2.0_dp, -2.0_dp, threshold) - ! Get next P1 -! call bml_multiply(B_bml, P0_bml, C_bml, 1.0_dp, 1.0_dp, threshold) -! call bml_multiply(Inv_bml(i), C_bml, P1_bml, 1.0_dp, 0.0_dp, threshold) + ! do i = 1, nsteps-1 + ! D = A^-1*P0 + ! call bml_multiply(Inv_bml(i), B0_bml, C0_bml, 1.0_dp, 0.0_dp, threshold) + ! call bml_multiply(C0_bml, B0_bml, B_bml, 1.0_dp, 0.0_dp, threshold) + ! B0 = A^-1*P0^2 + ! call bml_copy(B_bml,B0_bml) + ! B = I + D -P0*D + ! call bml_add(B_bml, C0_bml, -1.0_dp, 1.0_dp, threshold) + ! call bml_scale_add_identity(B_bml, 1.0_dp, 1.0_dp, threshold) + ! P1 = 2D*P1(I+D-P0*D) + ! call bml_multiply(C0_bml, P1_bml, C_bml, 1.0_dp, 0.0_dp, threshold) + ! call bml_multiply(C_bml, B_bml, P1_bml, 2.0_dp, 0.0_dp, threshold) + ! enddo + ! call bml_multiply(B0_bml, P1_bml, C_bml, 2.0_dp, 0.0_dp, threshold) + ! call bml_copy(P1_bml, B_bml) + ! call bml_add(B_bml, C_bml, 2.0_dp, -2.0_dp, threshold) + ! Get next P1 + ! call bml_multiply(B_bml, P0_bml, C_bml, 1.0_dp, 1.0_dp, threshold) + ! call bml_multiply(Inv_bml(i), C_bml, P1_bml, 1.0_dp, 0.0_dp, threshold) ! dPdmu = beta*P0(I-P0) @@ -1035,7 +1035,7 @@ subroutine prg_conjgrad(A_bml, p_bml, p2_bml, tmp_bml, d_bml, w_bml, cg_tol, thr do while (r_norm_new .gt. cg_tol) - ! write(*,*) r_norm_new + ! write(*,*) r_norm_new k = k + 1 if (k .eq. 1) then call bml_copy(tmp_bml, d_bml) @@ -1173,8 +1173,8 @@ subroutine prg_test_density_matrix(ham_bml, p_bml, beta, mu, nocc, osteps, occEr trdPdmu = beta * trdPdmu occErr = abs(trP0 - nocc) if (occErr .gt. occErrLimit) then - mu = mu + (nocc - trP0)/trdPdmu - end if + mu = mu + (nocc - trP0)/trdPdmu + end if !write(*,*) "mu = ", mu enddo diff --git a/src/prg_syrotation_mod.F90 b/src/prg_syrotation_mod.F90 index 44d7f91e..16f28b7e 100644 --- a/src/prg_syrotation_mod.F90 +++ b/src/prg_syrotation_mod.F90 @@ -190,7 +190,7 @@ subroutine prg_rotate(rot,r,verbose) v2(2)=rot%v2(2) v2(3)=rot%v2(3) - vQ(1)=rot%vQ(1) !Rotation center + vQ(1)=rot%vQ(1) !Rotation center vQ(2)=rot%vQ(2) vQ(3)=rot%vQ(3) @@ -232,7 +232,7 @@ subroutine prg_rotate(rot,r,verbose) v2=pq2-vQ endif - vtr(1)=0.0_dp !Translation + vtr(1)=0.0_dp !Translation vtr(2)=0.0_dp vtr(3)=0.0_dp diff --git a/src/prg_xlbokernel_mod.F90 b/src/prg_xlbokernel_mod.F90 index 81496915..8b36f9fb 100644 --- a/src/prg_xlbokernel_mod.F90 +++ b/src/prg_xlbokernel_mod.F90 @@ -332,8 +332,8 @@ subroutine prg_kernel_multirank(KRes,KK0_bml,Res,FelTol,L,LMAX,HO_bml,mu,beta,RX enddo !$OMP END PARALLEL DO - ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc,TIMERATIO,Max_Nr_Neigh) - ! Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k + ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc,TIMERATIO,Max_Nr_Neigh) + ! Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k row1 = 0.0_dp do J = 1,HDIM @@ -557,7 +557,7 @@ subroutine prg_full_kernel_latte(KK,DO_bml,mu0,RXYZ,Box,Hubbard_U, & enddo write(*,*) 'ewaldracc =', ewaldracc - write(*,*) 'ewaldkacc =', ewaldkacc + write(*,*) 'ewaldkacc =', ewaldkacc write(*,*) 'implcit response time =', respacc call Invert(JJ,KK,Nr_atoms) @@ -648,8 +648,8 @@ subroutine prg_full_kernel(KK,DO_bml,mu0,RX,RY,RZ,LBox,Hubbard_U,Element_Type, & Coulomb_Pot_Real(I) = Coulomb_Pot_Real_I enddo !$OMP END PARALLEL DO - ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc, & - ! TIMERATIO,Max_Nr_Neigh) + ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc, & + ! TIMERATIO,Max_Nr_Neigh) Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k diagonal = 0.0_dp @@ -766,8 +766,8 @@ subroutine prg_kernel_matrix_multirank(KRes,KK0_bml,Res,FelTol,L,LMAX,HO_bml,mu, enddo !$OMP END PARALLEL DO - ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc,TIMERATIO,Max_Nr_Neigh) - ! Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k + ! call Ewald_k_Space(Coulomb_Pot_k,Coulomb_Force_k,RX,RY,RZ,LBox,dq_v,Nr_atoms,Coulomb_acc,TIMERATIO,Max_Nr_Neigh) + ! Coulomb_Pot_dq_v = Coulomb_Pot_Real+Coulomb_Pot_k row1 = 0.0_dp do J = 1,HDIM diff --git a/tests/src/main.F90 b/tests/src/main.F90 index 105fa49a..9e437799 100644 --- a/tests/src/main.F90 +++ b/tests/src/main.F90 @@ -32,7 +32,7 @@ program main implicit none integer :: norb, mdim, verbose - type(bml_matrix_t) :: inv_bml(10) + type(bml_matrix_t) :: inv_bml(10) type(bml_matrix_t) :: ham_bml type(bml_matrix_t) :: rho_bml, rho1_bml type(bml_matrix_t) :: rho_ortho_bml @@ -215,7 +215,7 @@ program main error stop endif - case("prg_implicit_fermi_save_inverse") + case("prg_implicit_fermi_save_inverse") mu = 0.2_dp beta = 4.0_dp !nocc,osteps,occerrlimit @@ -224,7 +224,7 @@ program main do i = 1,norecs call bml_identity_matrix(bml_type,bml_element_real,dp,norb,norb,inv_bml(i)) - enddo + enddo call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,rho_bml) call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,ham_bml) @@ -232,10 +232,10 @@ program main call bml_read_matrix(ham_bml,'hamiltonian_ortho.mtx') call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,rho1_bml) - - mu = 0.2_dp + + mu = 0.2_dp call prg_implicit_fermi_save_inverse(inv_bml,ham_bml,rho_bml,norecs,nocc,mu,beta,1e-4_dp, threshold, 1e-5_dp, 1,occiter) - call prg_test_density_matrix(ham_bml,rho1_bml,beta,mu,nocc,1,1e-4_dp,threshold) + call prg_test_density_matrix(ham_bml,rho1_bml,beta,mu,nocc,1,1e-4_dp,threshold) write(*,*) mu call bml_scale(0.5_dp,rho_bml) call bml_add(rho1_bml,rho_bml,1.0_dp,-1.0_dp,threshold) @@ -1132,7 +1132,7 @@ program main call bml_add_deprecated(-1.0_dp,rho_bml,1.0_dp,rho1_bml,0.0_dp) error_calc = bml_fnorm(rho_bml) write(*,*)error_calc - + if(error_calc.gt.error_tol)then write(*,*) "Error is too high", error_calc