PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-02 06:51:37 +02:00
Code Releases Activity

Solved RESTART reading bug

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-02-06 22:41:35 +01:00
parent 1f8757afb6
commit 799d44324f
3 changed files with 26 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/HF/HFB.f90
View File

@ -144,7 +144,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
do iorb=1,nOrb do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * & P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo enddo
endif endif
@ -158,7 +158,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
do iorb=1,nOrb do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * & P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo enddo
endif endif
@ -343,11 +343,13 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
! Compute dipole moments, occupation numbers and || Anomalous density|| ! Compute dipole moments, occupation numbers and || Anomalous density||
! also print the restart file ! also print the restart file
eigVEC(:,:) = 0d0 deallocate(eigVEC,eigVAL)
allocate(eigVEC(nOrb,nOrb),eigVAL(nOrb))
eigVEC(1:nOrb,1:nOrb) = 0d0
eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb) eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb)
call diagonalize_matrix(nOrb2,eigVEC,eigVAL) call diagonalize_matrix(nOrb,eigVEC,eigVAL)
Occ(1:nOrb) = eigVAL(nOrb+1:nOrb2) Occ(1:nOrb) = eigVAL(1:nOrb)
c(1:nBas,1:nOrb) = matmul(X(1:nBas,1:nOrb),eigVEC(1:nOrb,nOrb+1:nOrb2)) c = matmul(X,eigVEC)
norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2))) norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2)))
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot) call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot)

10
src/HF/read_restart_HFB.f90
View File

@ -22,6 +22,7 @@ subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
integer :: nOrb_ integer :: nOrb_
double precision :: chem_pot_ double precision :: chem_pot_
double precision :: max_diff double precision :: max_diff
double precision :: val_read
double precision,allocatable :: eigVAL(:) double precision,allocatable :: eigVAL(:)
double precision,allocatable :: c_tmp(:,:) double precision,allocatable :: c_tmp(:,:)
double precision,allocatable :: S_mol(:,:) double precision,allocatable :: S_mol(:,:)
@ -37,11 +38,18 @@ subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb)) allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb))
c_tmp=0d0
S_mol=0d0
X_mol=0d0
open(unit=iunit,form='unformatted',file='hfb_bin',status='old') open(unit=iunit,form='unformatted',file='hfb_bin',status='old')
read(iunit) nBas_,nOrb_ read(iunit) nBas_,nOrb_
read(iunit) chem_pot_ read(iunit) chem_pot_
do iorb=1,nOrb do iorb=1,nOrb
read(iunit) c_tmp(:,iorb) do ibas=1,nBas
read(iunit) val_read
c_tmp(ibas,iorb) = val_read
enddo
enddo enddo
do iorb=1,nOrb do iorb=1,nOrb
read(iunit) eigVAL(iorb) read(iunit) eigVAL(iorb)

11
src/HF/write_restart_HFB.f90
View File

@ -19,6 +19,7 @@ subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
! Local variables ! Local variables
integer :: ibas,iorb,iunit=666 integer :: ibas,iorb,iunit=666
double precision :: val_print
! Dump results ! Dump results
@ -26,10 +27,16 @@ subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
write(iunit) nBas,nOrb write(iunit) nBas,nOrb
write(iunit) chem_pot write(iunit) chem_pot
do iorb=1,nOrb do iorb=1,nOrb
write(iunit) c(:,nOrb+1-iorb) do ibas=1,nBas
val_print = c(ibas,nOrb+1-iorb)
if(abs(val_print)<1d-8) val_print=0d0
write(iunit) val_print
enddo
enddo enddo
do iorb=1,nOrb do iorb=1,nOrb
write(iunit) Occ(nOrb+1-iorb) val_print = Occ(nOrb+1-iorb)
if(abs(val_print)<1d-8) val_print=0d0
write(iunit) val_print
enddo enddo
write(iunit) iunit write(iunit) iunit
close(iunit) close(iunit)