I am trying to calculate the electric field gradient (EFG) at the position of a certain nucleus for an excited electronic state. Unfortunately I do not find a way to do so, do you have any ideas?
I know how to calculate the EFG at the positions of the nuclei for the ground state using the $magnet section of the input file. One could also calculate the EFG via integration of the electron density and the multipole moment of interest (or integrating with 1/r and differentiating twice). However, I do not know if that is the way QChem calculates the EFG internally.
As far as the excited state calculation is concerned, currently the code reads in the last saved MO coefficients and density matrix from disk. If you know that you can get the right density saved to disk, you can run a job where you provide an existing scratch folder with the correct density from a previous job (it’s always stored in file 54) for a magnetman job with skip_scfman = true.
If you can give more detail about what kind of excited state calculations you’re doing, I can tell you whether or not this will work.
I am attempting to run a cis or cis(d) calculation on C6H6 + B (in the center of the molecule) and I am interested in the EFG of the excited states.
How do I save and read the MO coefficients of the excited states from disk, I have been trying to read in the density from the last job, but it always seems to be the ground state density. I ran
Would it also be possible to calculate the EFG according to equation 10.68 externally? But how would I obtain the MOs in the atomic orbital basis and the corresponding P^(a+b), so that I could integrate it myself?
Hi Johannes,
One quibble with your question is that “MO coefficients of excited states” is slightly ambiguous. There are amplitudes for a transition density but not true MOs in the excited state because the wave function is multiconfigurational. The closest thing are natural transition orbitals, which are eigenfunctions of the CIS/TDDFT 1-PDM. Those amplitudes must be in the formatted checkpoint file because IQmol reads them for excited-state plotting.
A partial answer to your question, or maybe an alternative solution, is that you can easily get the electrostatic potential for each excited state, printed out on a real-space grid in cube file format. Here’s a sample input for H2O:
$molecule
0 1
O
H 1 0.95
H 1 0.95 2 104.5
$end
$rem
method b3lyp
basis 6-31G*
esp_grid -3 ! excited state ESP
make_cube_files true ! generate cube files using plots grid
state_analysis true ! optional for this example, provides some extra analysis
cis_n_roots 2
cis_triplets false
$end
This will create a directory in the same place as your input/output files, with fairly obvious names for the various cube files that are created. There’s an option in the manual for electric field (ESP_EFIELD) but that doesn’t appear to have been ported to excited states.
You are right, please excuse my ambiguous question. I was just wondering how I could replicate the internal steps necessary to compute the EFG of excited states at the position of a certain nucleus. I still do not know how to obtain the relevant information (if that is even possible) to do so.
I have already tried the alternative solution you suggested. Unfortunately there seems to be a problem. As far as I know the electrostatic potential that is calculated via the esp_grid tag, is the total electrostatic potential including electron and nuclei terms. Since I am interested in the EFG at the position of a certain nucleus, I do not want this term to be included, however, I cannot simply subtract this ~1/r term of the potential of the nucleus since it is divergent at the position of interest, but the electric potential is not. I do not see how can work around this issue, but maybe you have something in mind?
Right, that makes sense. The existing code does have the ability to compute E-field at the nuclei (using ESP_EFIELD $rem variable). That is used by the Q-Chem/CHARMM interface for QM/MM gradient in the ground state. I took a look at the code and the electronic and nuclear contributions are computed separately in an obvious way, so could be printed separately without tremendous difficulty. That’s still not the E-field gradient but would that info be sufficient for finite-differencing? For E-field it’s ground-state only at this point (unlike ESP that has functionality for excited states), but that would be relatively easy to add.
Are you a developer? If so I could walk you through the code, send me an email. Otherwise you should submit as a feature request to Q-Chem. Not sure that I personally have the time to do it right now but it should be fairly straightforward, I think.
It looks like some bits are there, but require some hacking to bring together.
Run the CIS or CIS(D) calculation with state_analysis = true and gui = 2 in order to get the state densities in the formatted checkpoint file.
Run a ground state EFG calculation with dump_ao_integrals = true in the $magnet section.
The nuclear contribution to the EFG should be the same for an excited state as the ground state, since it’s the 2nd derivative of $V_{NN}$ divided by the nuclear charge.
A file will be produced called integrals.AO_EFG_1e.dat which is exactly what it says on the tin. It’s in Armadilloarma_ascii format, which is easy to parse with Python. For example, when adding this to the magnet_glycine_rhf sample job, it produces a file with the header
ARMA_CUB_TXT_FN008
80 80 60
which indicates that it’s a rank three object (a cube), where the first two dimensions are the number of basis functions, and the last dimension is a compound index (i*6) + j, where i is over all nuclei and j is over unique Cartesian components: XX, XY, XZ, YY, YZ, ZZ. (The matrices are symmetric with respect to Cartesian component interchange.) Glycine has 10 atoms, so 6*10 gives the total length along the last dimension. The file layout is such that each 80-by-80 matrix appears one after the other.
The electronic contribution to the EFG is going to be the contraction (elementwise multiplication and then complete accumulation) of the state density with each EFG matrix. For a restricted reference, you may need to multiply by two.
The total EFG tensor is the sum of the electronic and nuclear components. Remember that the tensor is symmetric.
I know it has been a while since this topic was discussed. I now have implemented the suggested approach, and it works in principle, but I obtain some strange results that are not consistent with the internal calculation of the EFG tensor in qchem.
I obtain the (electronic) EFG tensor in two ways: once from qchem directly by setting electric=true in the magnet section, and once by following the suggested steps. But the results differ. For small molecules or atoms it seems like it could be caused by numerics (since there are errors of 1e-4 at components that should be exactly zero) but for a larger molecule I obtain completely different results for some atoms.
I checked that the density is correctly normalized also considering the overlap matrix. Also, the summation is performed as suggested. So I am not sure what could cause these deviations, as I think that internally the EFG is calculated exactly the same way and the density and integrals are saved with 9 digit accuracy. Do you have any idea how I could find the cause for this problem?