Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Code review hack! #1

Draft
wants to merge 75 commits into
base: source_master
Choose a base branch
from
Draft

Code review hack! #1

wants to merge 75 commits into from

Conversation

ajwt
Copy link
Owner

@ajwt ajwt commented Feb 12, 2022

No description provided.

… UHF subroutines (now just duplicates of existing subroutines) for future development
…limiting steps are due to asymmetrised integrals, need to make slices of them next.
…specification of number of error matrices, documentation
Copy link
Owner Author

@ajwt ajwt left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Generally very nicely written code! I'd be interested in comparisons of performance between different scheduling or vectorization regimes.


file (GLOB_RECURSE sources src/*.f90 src/*.F90)
add_executable(els.x ${sources})
target_link_libraries (els.x /opt/OpenBLAS/lib/libopenblas.a -lpthread)
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't hurt to make this more configurable - having to have a dirty CMakeLists.txt hanging around because you don't have access to /opt isn't great practice.

src/hf.f90 Outdated
! ovlp%W holds the eigenvalues s_i of S
associate(diag=>ovlp%W)
do i = 1, size(ortmat, dim=1)
ortmat(i,i) = 1/sqrt(diag(i))
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For large bases this might end up blowing up owing to near linear dependence. Might be worth at least checking the value of diag(i) isn't too small, and throwing an error if it is.

src/hf.f90 Outdated
end do

if (.not. conv) then
write(iunit, '(1X, A)') 'Convergence not reached, please increase maxiter.'
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this worth throwing an error or perhaps returning something to note you haven't converged?

src/hf.f90 Outdated
end if

call deallocate_diis_t(diis)
deallocate(st%density_old)
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doing this here seems a bit out of place. Perhaps have a routine in system.F which does the allocation and deallocation of the bits of state_t?

src/hf.f90 Outdated
integer :: n_errmat = 6
integer :: n_active = 0
integer :: iter = 0
real(p), allocatable :: e(:,:,:)
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

definitely document what these are - also helpful to loosely know what the order of the subscripts is eg:
!F(norb,norb,n_max_diis_steps)

src/ccsd.f90 Outdated
end subroutine update_cc_energy

subroutine do_ccsd_t_spinorb(sys, int_store)
! Spinorbital formuation of CCSD(T)
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reference for the maths you're implementing?

@brianz98
Copy link

@ajwt I've made some significant changes to the code in the past 5 days, namely I've added the following functionalities:

  • Spin-free CCSD (by Piecuch et al.)
  • CCSD[T]/(T)
  • Renormalised CCSD[T]/(T)
  • Completely renormalised CCSD[T]/(T)
  • (these are all from the same paper, so I thought it'd be nice to code up and test all of them)
  • A python 'wrapper' around the code (basically just generates integrals from Psi4 and then runs it and captures the output and analyses the output, all in one script).

A head-scratcher for me is that the spin-free CCSD gives the wrong energies (and consequently I can't test the (T)/[T] methods). I've been at it for a week now, and completely coded the equations up from scratch again, and it's still the same, wrong energy. The extent of the error is about 0.5 mHa at eqm geometry for water/def2-svp, and almost 10 mHa at even moderately stretched 1.6 A, before completely diverging at 1.8 A. I've ruled out DIIS being the culprit as turning it off doesn't help.

In my repo, under src/ccsd.f90 I've added two new subroutines update_restricted_intermediates_debug (line 1216) and update_amplitudes_restricted_debug (line 1364), and a debug switch in do_ccsd_spatial (line 267). These subroutines are meant to be a no-frills, reference implementation of the equations (eq. 43-44 for the amplitude update equations, and table 1 for the intermediates) in Piecuch et al., as close to the original form as possible, with swapping of index order only when an MO slice is unavailable, using only basic do loops without any parallelisation. With these two reference subroutines I obtain the same wrong energy (and reassuringly slower). Could you please take a look at just those two subroutines, plus do_ccsd_spatial and init_cc, as I really need a fresh pair of eyes. Thanks so much!

p.s. the water/def2-svp data for r_OH = {1.00..0.20..1.80} can also be found in my repo under dat/, each with a reference.dat that is the correct energies from Psi4.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants