SCF Vibrational Circular Dichroism (VCD) Simulations#95
Conversation
|
Travis builds failing aren't necessarily a red flag. If the failures are caused by CEPA tests, it's because my tests require a newer version of Psi than what P4N is using. Lori tells me RHF_symmetry failures have been fixed in a local version of Psi, so that's also not your fault. |
|
@JonathonMisiewicz @loriab Can we mark the failing tests as xfail for now until Psi is updated? |
|
#96 now marks those as xfail. |
|
after rebase, a clean test suite is within your control. |
|
This is looking fairly close to mergeable. Do you have a list of items that you wish to accomplish here? |
|
@dgasmith Working on matching the final rotatory strengths now. One more commit and it should be ready for review. |
|
Here's where I'm at. I'm able to compute the APT's and match them perfectly with DALTON's values. However, my AAT's only match to about 3 decimal places when compared to DALTON's. As a result, when you contract the APTs and AATs together to get the rotatory strengths, the error is propogated and I match their values with even less precision. This seems problematic to me. I've broken it down, and I can match the nuclear contribution to the AAT's exactly. So the problem most definitely lies in the electronic contribution to the AAT's. I've broken that down even further. There are three quantities that make up the electronic contribution to the AAT's:
I know there are other codes that can compute VCD at the SCF level, viz. Gaussian. However, I can't figure out how to turn off GIAOs in Gaussian and therefore, until GIAOs are available in Psi4, this comparison is useless. As a further test, I've carried out the same VCD calculation with GIAOs turned on for both Gaussian and DALTON and their AATs don't match perfectly either. My natural inclination is that there must be a bug in my code somewhere, but seeing that their values don't match perfectly makes me wonder. At this point, I'm at a loss for what to do next to further nail down DALTON's AAT's and rotatory strengths. Any ideas? For comparison, my final rotatory strengths are: [50.53768056395719 -53.527791951020774 28.606794550012008 -14.649888335885528 1.0984924678221817 101.37800368356949] NOTE: I realize that my rotatory strengths have the opposite sign of DALTON's. This is also worrisome. However, when I did my tests with Gaussian vs. DALTON, Gaussian also had the same sign discrepancy when compared to DALTON's rotatory strengths. In other words, my rotatory strengths have the same sign as Gaussian, but both Gaussian and my values have opposite sign when compared to DALTON. TDC has brought this up with some higher ups at DALTON because we're fairly certain it is a problem on their end. Since our APT's and AAT's match signs, it doesn't make sense that an arbitrary negative sign would should up in each rotatory strength when they're simply contracted together. |
|
On another note, I don't know if anyone is still working on #68, but I've calculated IR intensities on the way in this code, and I'd be willing to clean that portion of the code up and use it as a reference implementation for that as well. |
|
This one is outside of my experience unfortunately. Pinging a few people to see if they have bandwidth to help out. PS: IR intensities would be awesome! Edit: Just to double check, everything is canonical integrals, Schwarz thresholds are off, integral thresholds are at max, basis set has been checked to be identical, SCF convergence is very high, etc? |
|
Great work Kirk! 👍 There could be an issue with the orbital connection for the response equations in DALTON. You should try forcing the use of the symmetric, rather than the natural connection: To make a reasonable comparison with Gaussian, try to look at convergence with the size of the basis set, as you push it to something stupidly large: the calculations with and without London orbitals should converge to the same results then. I am not sure I understand the sign issue, I'll think more on it and also review the actual code within tomorrow. |
|
Thanks @robertodr! |
@dgasmith Okay, great. How should I go about that? How do I push changes from my branch to a previously created PR?
Nuclear repulsion energies match identically to 15 decimal places, SCF energy matches to 8 decimal places. |
@robertodr Unforunately, I won't be able to crank up the basis set on my Python code very high since it's not very optimized and keeps everything in memory |
|
@kcpearce For IR intensities you would need to make a PR to the fork's branch. A bit complex, I would just pull that PR's branch to your fork and make a new PR. git remote add prbranch https://github.com/dsirianni/psi4numpy.git
git fetch prbranch IRintensities:newbranch
git checkout newbranch
git push origin newbranch |
|
Roberto has probably gotten everything technical already, but I'd still like to run this, and neither Psi4 1.3.2 nor 1.4a2.dev191 are running, getting so my two high-level comments so far are:
|
|
1.4a2.dev191 is very old. If you dont get an update offered from conda, try requesting the latest version explicitly. |
|
Here's an update on my testing of my code. I mentioned previously that the electronic component of my AAT was not matching DALTON's result exactly. So I decided to compare my results to an old paper by Roger Amos (https://www.sciencedirect.com/science/article/abs/pii/0009261487800465) in which he reports SCF APT's and AAT's using a 4-31G basis set for ehtylene oxide without the use of GIAO's. It was somewhat difficult to recreate the exact geometry that Amos used, but I feel like I've gotten rather close. Below are the results for my code, DALTON's code, and Amos's reported values for the AAT in units of T^-1 Angstrom^-1. Considering the fact that the geometry is most certainly slightly different than the one Amos uses, and I have no information as to what convergence criteria he was using, I feel rather good that I match Amos's values within +/-0.003 (with the exception of one single value corresponding to H1_y B_x - I'm very confused by this one). Unfortunately, I was hoping this would help me figure out why my code and DALTON's were differing slightly, but I still have no insight into why this might be the case. `MY CODE ELEC AAT (T^-1 A^-1): DALTONS ELEC AAT (T^-1 A^-1): AMOS ELEC AAT (T^-1 A^-1): ` |
|
Regarding the single H1_y B_x value that I was comparing to Amos's paper and getting a larger than expected error: I think this could be a simple fudging of numbers in his paper. My value for that element is That being said, if we assume that he did in fact have a typo for that single value, then I'm very happy with my comparisons to all of the values in Amos's paper. I know that my AATs differ slightly from DALTONs, and I'm going to continue to look into the source of this difference. However, I'm not sure if we should move forward and merge this reference implementation for the time being or wait until there is a better explanation for the difference. Thoughts? |
|
It looks like excellent detective work and a very plausible match to me. |
|
I generally agree, would like @robertodr and @berquist to comment before merging however. |
|
Unfortunately, due to everything described above, it's going to be tough to implement a test for this code. @dgasmith What do you recommend in the way of that? |
|
UPDATE: I've been able to track down the sources of error between my code's results and DALTON's results. 1.) The basis sets differ slightly. DALTON is slightly less precise with their basis set coefficients than Psi4/BSE. In order to get a direct comparison, I added 2.) The convergence thresholds for the CPHF solvers in DALTON had default values of After making these two changes, I was able to nail down DALTON's AAT values. For example, DALTON's AATs for H2O2 with a STO-3G basis set are the following: and my code's AATs for the same molecule and basis set are as follows: DALTON's corresponding rotatory strengths for the aforementioned example are: and my rotatory strengths are: I know that my rotatory strengths and DALTONs still differ in their sign, but @robertodr has informed me that it is a known issue that DALTONs VCD rotatory strengths differ when compared to Gaussian. As such, I am comfortable moving forward knowing that the sign of my values differ from theirs. As an additional comparison, I also compared the APTs and AATs between my code and DALTON's using the ethylene oxide example discussed above. and my APTs in a.u. are Lastly, DALTON's AATs in a.u. for this example are and my AATs in a.u. are Since ethylene oxide is non-chiral, the rotatory strengths are trivially zero and aren't worth comparing. Considering all of this, I am very happy with my comparisons. I'm planning on cleaning up my code and making one last commit. Once that's in, I should finally be ready for a review + merge. |
|
@kcpearce really great detective work here! |
|
Finally ready for review + merge! (Let me know if you'd like me to add a test for this reference implementation) |
dgasmith
left a comment
There was a problem hiding this comment.
LGTM, please do add tests.
Lets wait on a review from @robertodr before merging.
There was a problem hiding this comment.
- The reference output from DALTON needs to be updated (for the cranked-up convergence thresholds)
- Is the cooked-up basis part of Psi4? I suggest you download the proper STO-3G from the basis set exchange in DALTON format and re-run DALTON with the proper basis. Section 28.4 of the manual explains how to use custom bases (or I can help you run the calculation) - As for the test, I think you can copy the tensors of interest from the DALTON output (frequencies, IR intensities, AATs, APTs, VCD strengths), paste them in your Python source, and compare with your results.
|
I'm failing TravisCI with this error: It appears as though the Psi4 version in Travis is not recent enough to have my changes for exporting the electric dipole derivatives. |
robertodr
left a comment
There was a problem hiding this comment.
LGTM! I don't know how to solve the CI failures though 😞
|
@dgasmith Any idea about how to update the version of Psi4 for the CI? |
|
@dgasmith It looks like the last two commits aren't triggering CI. Any idea why? |
|
We probably ran out of CI bucks. Everything that is still on Travis CI OSS should be moved to GHA or similar since their policies changed. |
|
Finally passing! Everything is good to go on my end. |
|
Cool, thanks again for such an awesome addition! |
|
Hi All, It is not a good idea to compare with Dalton, because Dalton gives wrong result. Let me show you an example where Dalton gives the wrong sign of the Cotton effect. Input mol file (prop.mol): BASIS 6-311G prop Generated by Open Babel. Check overall charge below. AtomTypes=5 Charge=0 NoSymmetry Angstrom Charge=6.0 Atoms=2 C -1.8848000000 2.3694000000 -0.7438000000 C -1.7787000000 2.1817000000 0.7076000000 Charge=8.0 Atoms=1 O -1.1103000000 1.2941000000 -0.1970000000 Charge=1.0 Atoms=1 H -2.8053000000 2.0387000000 -1.2239000000 Charge=6.0 Atoms=1 C -1.1206000000 3.4432000000 -1.4696000000 Charge=1.0 Atoms=5 H -2.6064000000 1.7464000000 1.2624000000 H -1.1269000000 2.8385000000 1.2798000000 H -0.2208000000 3.7137000000 -0.9130000000 H -0.8168000000 3.0954000000 -2.4607000000 H -1.7367000000 4.3383000000 -1.5973000000 Input dal file (vcd.dal): **DALTON INPUT .RUN PROPERTIES **WAVE FUNCTIONS .HF **PROPERTIES .VIBANA .NOCMC .VCD **END OF DALTON INPUT VCD calculation: dalton -mol prop.mol -dal vcd.dal -np 12 And here is the result: mode frequency dip. str. rot. str. -------------------------------------------------- 1 3211.43 50.237 5.001 2 3151.22 63.941 -9.816 3 3131.96 12.242 5.948 4 3123.44 49.374 9.053 5 3120.03 22.720 -7.254 6 3067.38 36.849 -1.679 7 1694.15 7.012 -3.460 8 1655.13 17.999 3.665 9 1642.79 14.303 -2.897 10 1600.44 34.560 -9.039 11 1572.19 14.950 -8.679 12 1407.83 17.521 2.784 13 1332.80 5.762 6.925 14 1321.35 2.768 1.915 15 1293.42 19.217 12.609 16 1251.26 31.259 9.218 17 1174.12 19.545 -4.381 18 1080.39 47.470 37.988 19 1035.22 14.964 -27.033 20 910.46 250.409 -18.446 21 859.45 44.023 -5.530 22 465.33 30.281 -2.281 23 414.54 102.604 24.688 24 311.20 11.119 -5.066 Compared to the result calculated with other QM software at the same theory level, the sign of the rotational strength is wrong! |
Description
Adds SCF VCD capabilities to Psi4NumPy
TODOs
Notable points (developer or user-interest) that this PR has or will accomplish.
Any questions for the community?
Status