search

sharc-md / sharc

The SHARC molecular dynamics (MD) program suite is an ab initio MD software package developed to study the excited-state dynamics of molecules.
https://www.sharc-md.org
GNU General Public License v3.0
64 stars 35 forks source link

Assistance Required for Simulated Transient Absorption Spectra for Specific Delay Times #135

Open lukhman9020 opened 1 month ago

lukhman9020 commented 1 month ago

Dear Sir, I am trying to simulate transient absorption spectra for some specific delay times as you did in your literature titled "Ultrafast Intersystem Crossing Dynamics of 6-Selenoguanine in Water." However, I have been unable to obtain the correct spectra. I can generate TAS spectra for energy vs. time (fs) as shown in the tutorial (2.9.10 Transient Spectra), but I need guidance on how to generate TAS for specific delay times. Sir , Could you please assist me in generating TAS for some specific delay times?

Best regards, LUKHMANUL HAKEEM K.

lukhman9020 commented 1 month ago

transient absorption spectra for some specific delay times

Dear Sir, Please consider my request and assist me in simulating transient absorption spectra for specific delay times. I look forward to your positive response.

Best regards, LUKHMANUL HAKEEM K.

maisebastian commented 1 month ago

Dear lukhman9020, sorry, I am currently traveling and have less time. So just briefly - in the past, for such TAS as you describe, we put together some quick Python plus Bash script solutions. These are not generally applicable and are not part of the SHARC package. Basically, you would have to: 1) Decide on the times you want to compute TAS for 2) For each traj and each time, make a folder 3) For each traj and time, get the coordinates from output.xyz, put them into the made folder 4) Add the necessary keywords and files in each folder to make a SHARC interface single point calculation. Here it is important that in your method, the number of states does not affect the results, because in the new calculation, the active state needs to be assignable. So TDDFT and ADC(2) would work, but CASSCF/CASPT2 might not. 5) Run all the single point calculations 6) For each traj and time, compute the energy differences and oscillator strengths from the active state to all other states. 7) Convolute the dE and fosc information to get a spectrum for each time step 8) Optionally, subtract the ground state spectrum if you are interested in ground state bleach.

Best, Sebastian

lukhman9020 commented 1 month ago

Dear lukhman9020, sorry, I am currently traveling and have less time. So just briefly - in the past, for such TAS as you describe, we put together some quick Python plus Bash script solutions. These are not generally applicable and are not part of the SHARC package. Basically, you would have to:

  1. Decide on the times you want to compute TAS for
  2. For each traj and each time, make a folder
  3. For each traj and time, get the coordinates from output.xyz, put them into the made folder
  4. Add the necessary keywords and files in each folder to make a SHARC interface single point calculation. Here it is important that in your method, the number of states does not affect the results, because in the new calculation, the active state needs to be assignable. So TDDFT and ADC(2) would work, but CASSCF/CASPT2 might not.
  5. Run all the single point calculations
  6. For each traj and time, compute the energy differences and oscillator strengths from the active state to all other states.
  7. Convolute the dE and fosc information to get a spectrum for each time step
  8. Optionally, subtract the ground state spectrum if you are interested in ground state bleach.

Best, Sebastian

Dear Sir,

Thank you very much for your valuable reply despite your busy schedule. I apologize for any inconvenience caused and wish you a happy and safe journey.

Sincerely, Lukhmanul Hakeem K

lukhman9020 commented 1 month ago

Dear lukhman9020, sorry, I am currently traveling and have less time. So just briefly - in the past, for such TAS as you describe, we put together some quick Python plus Bash script solutions. These are not generally applicable and are not part of the SHARC package. Basically, you would have to:

  1. Decide on the times you want to compute TAS for
  2. For each traj and each time, make a folder
  3. For each traj and time, get the coordinates from output.xyz, put them into the made folder
  4. Add the necessary keywords and files in each folder to make a SHARC interface single point calculation. Here it is important that in your method, the number of states does not affect the results, because in the new calculation, the active state needs to be assignable. So TDDFT and ADC(2) would work, but CASSCF/CASPT2 might not.
  5. Run all the single point calculations
  6. For each traj and time, compute the energy differences and oscillator strengths from the active state to all other states.
  7. Convolute the dE and fosc information to get a spectrum for each time step
  8. Optionally, subtract the ground state spectrum if you are interested in ground state bleach.

Best, Sebastian

Dear Sir,

I have successfully completed steps 1, 2, 3, 4, and 5. However, after completing the single-point calculations, I am confused about how to proceed with steps 6 and onward. Could you please explain which script I should use to complete these steps, and provide a detailed explanation of the process? I look forward to your positive response.

Best regards, LUKHMANUL HAKEEM K.

maisebastian commented 1 month ago

As I said before, computing TAS in the way you want to is not fully implemented in SHARC 3, so there is no script that you can directly use. Some brief hints:

lukhman9020 commented 1 month ago

As I said before, computing TAS in the way you want to is not fully implemented in SHARC 3, so there is no script that you can directly use. Some brief hints:

  • You can use $SHARC/QMout_print.py in Step 6. You would still need to identify the currently active state, in the indexing of the new single point calculation. E.g., if your trajectories were running with "3 0 3" states and the current MCH active state is 4 (T1), and if your single point was run with "10 0 10", then you would need to run QMout_print.py -s "10 0 10" -S 11, because the T1 is state 11 in the single point calculation.
  • For step 7, you will need to write a script that takes the outputs of QMout_print.py and convolute them. I suggest that you reuse the spectrum class from spectrum.py for this purpose.
  • Step 8 just involves subtracting two spectra. Make sure that they are both using the same settings for the energy axis (Emin, Emax, number of points).

Dear Sir, Thank you very much for your reply .

Best regards, LUKHMANUL HAKEEM K.

lukhman9020 commented 4 weeks ago

6. For each traj and time, compute the energy differences and oscillator strengths from the active state to all other states.

Dear Sir, After completing the single-point calculation in the SHARC-ORCA interface, I used the command "$SHARC/QMout_print.py -s "8 0 8" -S 9 " to compute the energy differences and oscillator strengths from the active state to all other states. However, I encountered an error as shown below: " Traceback (most recent call last): File "/home/mhgroup1/kumar/sharc/bin/QMout_print.py", line 443, in main() File "/home/mhgroup1/kumar/sharc/bin/QMout_print.py", line 301, in main qmoutfile = args[0] IndexError: list index out of range " Please help me to solve the error.

Best regards, LUKHMANUL HAKEEM K.

maisebastian commented 4 weeks ago

Sorry, you also need to give the file name of the QM.out as argument, like so: $SHARC/QMout_print.py -s "8 0 8" -S 9 QM.out

lukhman9020 commented 4 weeks ago

Sorry, you also need to give the file name of the QM.out as argument, like so: $SHARC/QMout_print.py -s "8 0 8" -S 9 QM.out

Dear Sir,

Thank you for your reply. I have one doubt, though it is not related to the above response. This is regarding how to extract the SOC value from the LVC.template file in cm^-1. When I inspect the LVC.template file, I can see the SOC R and SOC I matrices. How do we get the SOC value from this, for example, the SOC value between S1 and T1?

Here, I have parametrized the LVC using the SHARC-ORCA interface. Should I inspect the orca.log file in the scratch directory of DSPL_000_eq to get the SOC value?

Best regards, LUKHMANUL HAKEEM K.

maisebastian commented 4 weeks ago

This depends on what kind of SOC value you want to obtain.

If you are interested in the SOC values of the diabatic states, then you can take them from the LVC.template file. There, the SOC matrix is split into real and imaginary parts, and the full matrix for all MS components is present. So if you have 3 singlets and 2 triplets, then the SOC matrices are 3+3*2=9-dimensional. The order is S0, S1, S2, T1, T2, T1, T2, T1, T2. The SOC matrices are given in Hartree. To convert to cm-1, multiply by 219474. To get the effective SOC for S-T, compute the complex vector norm of the three MS components of the SOC <Sn|Hsoc|Tm_MS>.

If you are interested in the SOC values of the adiabatic states at a certain geometry, make a single point calculation with SHARC_LVC.py and get the Hamiltonian from the QM.out file. The ordering and conventions are the same as given above.

lukhman9020 commented 4 weeks ago

This depends on what kind of SOC value you want to obtain.

If you are interested in the SOC values of the diabatic states, then you can take them from the LVC.template file. There, the SOC matrix is split into real and imaginary parts, and the full matrix for all MS components is present. So if you have 3 singlets and 2 triplets, then the SOC matrices are 3+3*2=9-dimensional. The order is S0, S1, S2, T1, T2, T1, T2, T1, T2. The SOC matrices are given in Hartree. To convert to cm-1, multiply by 219474. To get the effective SOC for S-T, compute the complex vector norm of the three MS components of the SOC <Sn|Hsoc|Tm_MS>.

If you are interested in the SOC values of the adiabatic states at a certain geometry, make a single point calculation with SHARC_LVC.py and get the Hamiltonian from the QM.out file. The ordering and conventions are the same as given above.

Dear Sir, Thank you for your reply. Best regards, LUKHMANUL HAKEEM K.

lukhman9020 commented 3 weeks ago

Sorry, you also need to give the file name of the QM.out as argument, like so: $SHARC/QMout_print.py -s "8 0 8" -S 9 QM.out

Dear Sir,

I successfully extracted results from all QM.out files and stored them into a spectrum file. However, when I tried to run the command $SHARC/spectrum.py -o spectrum.out spectrum, I encountered the following error: " Number of grid points: 500 Energy range: 0.000 to 10.000 eV Lineshape: Gaussian (FWHM=0.100 eV) File malformatted! " Below is the content of the spectrum file: Number of states: [8, 0, 8] State Label E (E_h) dE (eV) f_osc Spin

1          S00 -6006.7486710000  -2.35386740  -0.00000000   1.0000
2          S01 -6006.6315550000   0.83302247   0.00000000   1.0000
3          S02 -6006.5983890000   1.73551563   0.00000000   1.0000
4          S03 -6006.5868630000   2.04915419   0.00000000   1.0000
5          S04 -6006.5779010000   2.29302272   0.00000000   1.0000
6          S05 -6006.5760800000   2.34257467   0.00000000   1.0000
7          S06 -6006.5680220000   2.56184410   0.00000000   1.0000
8          S07 -6006.5677060000   2.57044290   0.00000000   1.0000
9          T01 -6006.6621680000   0.00000000   0.00000000   3.0000 #initial state

10 T02 -6006.6434510000 0.50931570 0.00000000 3.0000 11 T03 -6006.6163630000 1.24641800 0.00000000 3.0000 12 T04 -6006.6021380000 1.63350011 0.00000000 3.0000 13 T05 -6006.5989000000 1.72161061 0.00000000 3.0000 14 T06 -6006.5893150000 1.98243184 0.00000000 3.0000 15 T07 -6006.5828670000 2.15789092 0.00000000 3.0000 16 T08 -6006.5808990000 2.21144295 0.00000000 3.0000 1 S00 -6001.3172260000 2.23313043 0.00000000 1.0000 2 S01 -6001.3375940000 1.67888872 0.00000000 1.0000 ......................................................................................................................................... .......................................................................................................................................... Could you please help me resolve the "File malformatted!" error?

Additionally, I am curious about generating a Python script to create the TAS spectrum. I need clarification on some information, specifically, image

what is the "E" in the equation? Is it similar to "E_h" in the spectrum file?. Please help me .

Best regards, LUKHMANUL HAKEEM K

maisebastian commented 3 weeks ago

The script spectrum.py is only made to read excited-state data from SHARC's initconds format. The output of QMout_print.py is not in that format, therefore it does not work.

You can use the attached script (extract the .py file) to convolute excitation energies and oscillator strengths, as is given by QMout_print.py. conv_energy.zip The command should look somewhat like this: python conv_energy.py -E 4 -F 5 -e -10 10 -N 1000 -f 0.5 -L g yourfile > spectrum This will simply read energies from column 4 and oscillator strengths from column 5 and convolute all the lines it can find. In this way, you can concatenate the results of all calculations at the same delay time into one long list of transitions, then convolute once for each delay time.

About your question on the equation, the "E" is simply the coordinate of the spectrum, which is a function of energy...

Best, Sebastian

lukhman9020 commented 3 weeks ago

The script spectrum.py is only made to read excited-state data from SHARC's initconds format. The output of QMout_print.py is not in that format, therefore it does not work.

You can use the attached script (extract the .py file) to convolute excitation energies and oscillator strengths, as is given by QMout_print.py. conv_energy.zip The command should look somewhat like this: python conv_energy.py -E 4 -F 5 -e -10 10 -N 1000 -f 0.5 -L g yourfile > spectrum This will simply read energies from column 4 and oscillator strengths from column 5 and convolute all the lines it can find. In this way, you can concatenate the results of all calculations at the same delay time into one long list of transitions, then convolute once for each delay time.

About your question on the equation, the "E" is simply the coordinate of the spectrum, which is a function of energy...

Best, Sebastian

Dear Sir,

Thank you for your prompt response and for sharing the script with me.

Best regards, LUKHMANUL HAKEEM K

lukhman9020 commented 3 weeks ago

python conv_energy.py -E 4 -F 5 -e -10 10 -N 1000 -f 0.5 -L g yourfile > spectrum

Dear Sir, I used the Python script provided to extract the dE and fosc values from my " init "file and also generated a file for plotting the absorption spectrum. However, the plot I obtained does not appear to be correct. Below, I have included the contents of the input file and the output file. Input File (init.excited): 1 S00 -6006.7486710000 -2.35386740 -0.00000000 1.0000 2 S01 -6006.6315550000 0.83302247 0.00000000 1.0000 3 S02 -6006.5983890000 1.73551563 0.00000000 1.0000 4 S03 -6006.5868630000 2.04915419 0.00000000 1.0000 5 S04 -6006.5779010000 2.29302272 0.00000000 1.0000 6 S05 -6006.5760800000 2.34257467 0.00000000 1.0000 7 S06 -6006.5680220000 2.56184410 0.00000000 1.0000 8 S07 -6006.5677060000 2.57044290 0.00000000 1.0000 9 T01 -6006.6621680000 0.00000000 0.00000000 3.0000 #initial state 10 T02 -6006.6434510000 0.50931570 0.00000000 3.0000 11 T03 -6006.6163630000 1.24641800 0.00000000 3.0000 12 T04 -6006.6021380000 1.63350011 0.00000000 3.0000 13 T05 -6006.5989000000 1.72161061 0.00000000 3.0000 14 T06 -6006.5893150000 1.98243184 0.00000000 3.0000 15 T07 -6006.5828670000 2.15789092 0.00000000 3.0000 16 T08 -6006.5808990000 2.21144295 0.00000000 3.0000 1 S00 -6001.3172260000 2.23313043 0.00000000 1.0000 2 S01 -6001.3375940000 1.67888872 0.00000000 1.0000 3 S02 -6001.3367610000 1.70155581 0.00000000 1.0000 4 S03 -6001.3360600000 1.72063100 0.00000000 1.0000 5 S04 -6001.3351710000 1.74482193 0.00000000 1.0000 6 S05 -6001.3236430000 2.05851491 0.00000000 1.0000 7 S06 -6001.3234600000 2.06349459 0.00000000 1.0000 8 S07 -6001.3216640000 2.11236626 0.00000000 1.0000 9 T01 -6001.3992920000 0.00000000 0.00000000 3.0000 #initial state 10 T02 -6001.3983240000 0.02634063 0.00000000 3.0000 11 T03 -6001.3951450000 0.11284566 0.00000000 3.0000 12 T04 -6001.3935090000 0.15736350 0.00000000 3.0000 13 T05 -6001.3727610000 0.72194555 0.00000000 3.0000 14 T06 -6001.3638700000 0.96388207 0.00000000 3.0000 15 T07 -6001.3629210000 0.98970569 0.00000000 3.0000 16 T08 -6001.3569530000 1.15210330 0.00000000 3.0000 ............................................................................................................................................

Output File (Used for Plotting the Spectrum):

0.000000000 nan 0.000000000 0.020000000 61950.000000000 0.000000000 0.040000000 30975.000000000 0.000000000 0.060000000 20650.000000000 0.000000000 0.080000000 15487.500000000 0.000000000 0.100000000 12390.000000000 0.000000000 0.120000000 10325.000000000 0.000000000 0.140000000 8850.000000000 0.000000000 0.160000000 7743.750000000 0.000000000 0.180000000 6883.333333333 0.000000000 0.200000000 6195.000000000 0.000000000 0.220000000 5631.818181818 0.000000000 0.240000000 5162.500000000 0.000000000 0.260000000 4765.384615385 0.000000000 0.280000000 4425.000000000 0.000000000 0.300000000 4130.000000000 0.000000000 0.320000000 3871.875000000 0.000000000 0.340000000 3644.117647059 0.000000000 0.360000000 3441.666666667 0.000000000 0.380000000 3260.526315789 0.000000000 0.400000000 3097.500000000 0.000000000 0.420000000 2950.000000000 0.000000000 0.440000000 2815.909090909 0.000000000 0.460000000 2693.478260870 0.000000000 0.480000000 2581.250000000 0.000000000 0.500000000 2478.000000000 0.000000000 0.520000000 2382.692307692 0.000000000 0.540000000 2294.444444444 0.000000000 0.560000000 2212.500000000 0.000000000 0.580000000 2136.206896552 0.000000000 0.600000000 2065.000000000 0.000000000 0.620000000 1998.387096774 0.000000000 0.640000000 1935.937500000 0.000000000 0.660000000 1877.272727273 0.000000000 0.680000000 1822.058823529 0.000000000 0.700000000 1770.000000000 0.000000000 0.720000000 1720.833333333 0.000000000 0.740000000 1674.324324324 0.000000000 0.760000000 1630.263157895 0.000000000 ..................................................................................................

image

Please help me .

Best regards, LUKHMANUL HAKEEM K

maisebastian commented 3 weeks ago

Two problems: 1) Your oscillator strengths are all zero in the output of QMout_print.py. I apologize that I did not realize earlier, but in regular TDDFT, one can usually not obtain the excited->excited oscillator strengths (this is possible in ADC(2)). In ORCA 5, there is an option to compute approximate excited->excited oscillator strengths (called "dotrans"). In https://github.com/sharc-md/sharc/blob/sharc3preview/bin/SHARC_ORCA.py, I have added this option. In order to use it, you have to request "dm" in QM.in and add "dotrans" to ORCA.template. The option is not compatible with SOCs and one gets only the excited->excited oscillator strengths for the singlet states. If that works for you, then you can use the mentioned script, rerun your single point calculations, and then continue as discussed above.

2) You need to plot columns 1 versus 3 or 2 versus 3. Column 1 is energy, column 2 is wave length (if the energy was given in eV), column 3 is the intensity.

Best, Sebastian

lukhman9020 commented 3 weeks ago

Dear Sir, Thank you for your reply. Best regards, LUKHMANUL HAKEEM K.

Dear Sir, Thank you for your reply. Best regards, LUKHMANUL HAKEEM K.

maisebastian commented 3 weeks ago

Dear user, unfortunately, the ORCA implementation of ADC(2) does not have gradients, as far as I know. Without gradients, running dynamics is not really possible. Therefore, so far we did not include ADC(2) in the options for the SHARC-ORCA interface. Hence, if you do not have access to Turbomole, you will only be able to use TDDFT as single-reference methods.

Best, Sebastian

lukhman9020 commented 3 weeks ago

Dear user, unfortunately, the ORCA implementation of ADC(2) does not have gradients, as far as I know. Without gradients, running dynamics is not really possible. Therefore, so far we did not include ADC(2) in the options for the SHARC-ORCA interface. Hence, if you do not have access to Turbomole, you will only be able to use TDDFT as single-reference methods.

Best, Sebastian

Dear Sir, Thank you for your reply. Best regards, LUKHMANUL HAKEEM K.